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PREFACE 


The  primary  purpose  of  this  report  is  to  provide  the  reader  with 
enough  background  in  scattering  theory  and  other  disciplines  to  enable 
him  to  better  understand  what  is  involved  in  calculations  of  the  energy 
eigenvalues  of  electrons  in  a  cluster  model.  This  report  discusses 
Hartree-Fock  equations,  self-consistent  field  theory,  and  scattering 
theory . 

The  author  supplemented  the  presentation  with  derivations  which  are 
ordinarily  left  out  of  journal  articles.  Although  not  an  easy  task,  the 
author  consolidated  material  from  many  references  in  a  moderately 
successful  manner.  The  author  wishes  to  acknowledge  all  of  the  authors 
cited  in  the  reference  section,  particularly  K.  H.  Johnson  whose  material 
Section  IX  was  based  upon. 


The  author  wishes  to  thank  Or.  Gust  Bambakidis  for  helping  him  to 
run  the  required  computer  programs.  Without  his  help,  the  computations 
could  not  have  been  made.  He  also  wishes  to  express  appreciation  to 
Dr.  Bambakidis  for  suggesting  this  report  topic  as  well  as  for  providing 
him  with  new  insights  into  a  complex  physical  problem.  The  author  wishes 
to  thank  Dr.  Donn  Shankland  of  the  Air  Force  Institute  of  Technology 
(AFIT).  Over  the  course  of  attending  his  lectures  on  Quantum  Mechanics 
and  Its  Formalism,  he  became  exposed  to  Green's  Functions  and  other 
concepts  which  provided  the  necessary  tools  for  the  theory  contained  in 
this  report. 


Finally,  the  author  wishes  to  express  his  sincerest  thanks  to  Mr. 
Frank  R.  Kile,  Jr.  for  providing  his  artistic  talents  to  this  project. 
Frank  drew  all  the  figures  in  the  text  in  the  most  professional  fashion. 
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SECTION  I 
INTRODUCTION 


This  report  describes  a  computer-based  calculation  of  the  diamond 
structured  complexes,  CgHg  and  C5H1 2  *  It  results  in  a  calculation  of 
the  energy  eigenvalues  of  electrons  in  these  two  molecular  complexes. 

The  computer  programs  are  based  on  the  self-consistent-field  scattered 
wave  +  X-ralpha  (SCF-SW+Xa)  method,  a  concept  originally  fathered  by  John 
C.  Slater  and  later  developed  by  Keith  Johnson  and  co-workers  at  M.I.T. 


The  report  is  divided  into  three  parts  (1)  theoretical  aspects, 
(2)  discussion  of  the  programs,  and  (3)  discussion  of  the  results. 


The  main  thrust  of  the  report  is  not  the  calculation  itself,  but  the 
development  of  the  theoretical  background  engendered  in  the  calculation. 
The  author  preferred  this  approach  because  the  theoretical  background  is 
not  readily  at  hand,  but,  on  the  contrary,  scattered  about  an  incoherent 
fashion  in  the  literature.  The  inception  of  X-a  scattered  wave  theory 
began  in  1937.  It  finally  reached  a  degree  of  maturity  in  the  late 
sixties.  Much  work  has  been  done  on  the  subject,  and  many  disciplines 
enter  into  the  fabric  of  the  topic:  Hartree-Fock  theory,  numerical 
methods,  quantum  mechanics,  quantum  statistical  mechanics,  theory  of 
solids,  group  theory,  and  mathematical  physics  are  all  interwined. 


It  is  the  author's  desire  that  this  report  provide  some  degree  of 
coherency  to  the  topic  as  well  as  insight  into  the  true  complexity  of 
the  problem. 
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SECTION  II 

INTUITIVE  DESCRIPTION  OF  THE  SELF-CONSISTENT  APPROACH 
TO  THE  ELECTRONIC  STRUCTURE  OF  ATOMS 

0.  R.  Hartree  devised  the  self-consistent-field  method  in  1928. 
Essentially,  he  assumed  that  each  electron  moves  in  an  averaged  potential 
arising  from  the  other  electrons  in  an  atomic  system.  The  Schrodinger 
equation  for  an  electron  moving  in  that  potential  is  then  solved;  this 
is  a  relatively  simple  problem  since  spherical  symmetry  is  assumed. 

The  wave  function  with  the  desired  quantum  number  in  that  potential 
is  then  selected.  This  wave  function  v  is  then  assumed  to  be  used  in 
finding  the  charge  density  arising  from  that  particular  electron  in 
question.  The  total  charge  density  arising  from  all  the  electrons 
present  in  the  atom  is  built  up.  The  potential  arising  from  this  charge 
density  (solution  of  Poisson's  equation)  is  calculated  by  electrostatics, 
and  then  the  requirement  for  self-consistency  is  made,  that  is,  this 
final  potential  must  agree  with  the  initial  one  which  had  been  assumed 
in  setting  up  Schrodinger' s  equation. 

Reiterating,  Hartree  found  that  he  could  set  up  a  very  manageable 
procedure  for  determining  the  self-consistent  field,  based  on  a  method 
of  iteration  or  successive  approximation.  He  started  by  assuming  trial 
wave  functions,  which  he  hoped  were  nearly  the  final  ones.  From  these 
functions  he  determined  charge  densities  and  potentials,  and  solved  the 
Schrodinger's  equation,  coming  out  with  final  wave  functions.  Since  he 
could  not  have  expected  to  have  achieved  self-consistency,  these  final 
wave  functions  would  not  agree  with  the  initial  ones.  But  then  he  started 
a  new  cycle  of  operation  using  the  final  functions  of  one  step  as  the 
starting  functions  of  the  next  step,  finding  charge  densities,  potentials, 
and  so  on,  from  them  and  proceeding  just  as  before.  Fortunately,  he 
found  that  this  process  converged  by  carrying  through  a  number  of  cycles; 
he  came  out  with  the  wave  functions  which  were  self-consistent  to  a  very 
good  approximation. 
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SECTION  III 

THE  HARTREE  SELF-CONSISTENT- FIELD  METHOD 

In  the  Hartree  Method,  it  is  assumed  that  each  electron  moves  in  a 
central  field,  produced  by  a  nucleus  and  the  spherically-averaged  potential 
fields  of  every  other  electron. 


We  assume,  first  of  all,  that  the  one-electron  wave  functions  or 
orbitals,  as  we  shall  call  them,  have  the  form 


...«  -  <-,)(v""‘l)/2  mZEEIK 


m„ 


X  RnH(r)  P£  1  <cos0)  expdm^), 
where  the  substitution 

Pnlt(r)=rfWr) 


will  be  made,  with  RnJl  or  P^  normalized  so  that 


|^[pnJl(r)]2dr-i 


The  function  Ps.  (-cose)  is  the  associated  Legendre  polynomial. 


Next,  we  inquire  how  a  wave  function  for  the  N-electron  atom  can 
be  constructed  out  of  these  one-electron  orbitals.  Each  electron  is 
supposed  to  move  quite  independently  of  the  others,  being  acted  on  by 
the  others  only  in  an  averaged  manner.  This  implies  that  the  wave 
function  \ i>  should  be  a  product  of  functions  of  the  various  electrons  or: 


-<J>  (>%e,j5...rN0N0N)  = 


Un1£1m1(r,0’0)*--Un 


N*Nmi 


(rn-VV 


M 


(1) 


abiil 


v.yv,v, 
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where  the  U's  on  the  RHS  of  Equation  1  are  the  orbitals,  each  one  a 
function  of  the  coordinates  of  a  single  electron.  We  assign  the  quantum 
numbers,  n.n.rn^,  etc.  for  the  various  electrons  according  to  the 
configurations  in  which  we  expect  the  atom  to  be  found. 

Equation  1  does  not  represent  an  exact  solution  of  Schrodinger 's 

equation  for  the  atom;  it  is  an  approximate  function  containing  the 

undetermined  radial  functions  R  (r)  as  quantities  may  be  varied.  For 

nNstN 

this  reason,  we  can  carry  through  this  as  a  variational  problem,  and 
the  wave  functions  prove  to  be  very  nearly  those  of  Hartree's  method. 

If  we  slightly  modify  the  calculation,  we  get  Hartree's  equation.  We 
then  modify  it  by  not  exactly  computing  the  average  energy,  but  by  making 
the  spherical  average  of  the  potential  of  an  electron  which  underlies 
Hartree's  approximation. 

1.  EXPECTATION  VALUE  OF  THE  ATOMIC  HAMILTONIAN 

The  first  step  in  applying  the  variation  method  to  Hartree's  wave 
function  is  to  compute  the  average  value  of  the  Hamiltonian  for  the  wave 
function  of  Equation  1  making  the  approximation  of  spherical  averaging. 
The  Hamiltonian  for  an  N-electron  system  is 

(H)op  *  -Z(i)V?  -  Z(1)  E(pairs  i,j)  -2- 

H  1  rij  (2) 

where  the  summation  over  i  is  also  over  all  N  electrons,  that  summation 
over  pairs  is  over  each  pair  counted  once,  excluding  the  case  i=j.  The 
first  term  in  Equation  2  is  the  sum  of  the  kinetic  energies  of  all 
electrons,  the  second  is  the  potential  energy  in  the  field  of  the 
nucleus,  r^  being  the  distance  from  the  nucleus  to  the  ith  electron, 
and  the  last  is  the  repulsive  Coulomb  potential  of  interaction  between 
the  ith  and  jth  electrons.  We  understand  that  we  are  to  perform  a 
spherical  average  of  this  repulsive  interaction.  Equation  2  is  in 
atomic  units.  It  is  convenient  to  rewrite  Equation  2  as: 

(h)qq  s  £(1)f.j  +  £(pairs  i ,j )  gM 

P  1  (3) 
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where 


fi  ="Vi  “  K"’  9ij 


We  now  let  the  operator  of  Equations  2  or  3  operate  on  the  wave 
function,  multiply  by  the  conjugate  of  the  wave  function,  and  integrate 
over  all  values  of  the  coordinates. 


First,  we  will  look  at  the  operator  f^ 

Vi.  <*,  ■  (df|1) 


iii. 


i  i  A,. 


Next,  we  will  consider  one  of  the  two-electron  operations  g...  It 

*  J 

operates  on  two  orbitals  with  indices  i  and  j.  We  can  then  integrate 

over  the  coordinates  of  all  electrons  except  these  two  and  find  as  the 

contribution  of  this  operation  g. .  to  the  average  Hamiltonian 

*  3 

JltViV  9ij  Vl^m  u„  i.m.  ]dvi  dvS  .  , 

*1  J  J  ^  1  1  J  J  A.Javerage 


=  (iJ|g|ij)av 

We  can  now  combine  these  results,  and  find  that 

(H)av  =  E(i )  (i | f | i )  +  E(pairs  i.j)  (ij|g|ij)  av 
2.  ENERGY  INTEGRALS  FOR  THE  HARTREE  CALCULATION 


We  know 

2  ,  (vH|)/2 

-Vu  * 


[gEmERDI 

(A+|mJ)l 


P£  ^  (cosO)  exp(1m  0) 
x  i  r  L  sJ- 


*  «  *  •  1 1 
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where  u.n.t.M^,  r,  e,  |  stand  for  ,  i . ,  etc.  if  we  are  dealing  with 
the  function  u,. .  We  multiply  by  UT  and  integrate  over  coordinates, 
making  use  of  the  properties  of  the  associated  Legendre  functions.  We 
f  i  nd : 


IViV  <‘V*)Vi'Y  dVl  "  lo  A,lf(ri)  [' 


4* 

dr1 


^(4+1 


-T 


5 


X  pn.Ai  dr.. 


(8) 


As  for  the  term  -  —  in  the  operator  f .. ,  we  can  carry  out  immediately 

the  integration  over  angles  which  leaves  a  radial  integration  only. 
Hence,  we  arrive  at 


Vi(r’)dr< 


(9) 


We  shall  denote  this  one-electron  integral  as 

(i  | f 1 1 )  =  Kn^^V  (10) 

Next,  we  consider  the  two-electron  integral  of  Equation  5.  This  represents 
the  Coulomb  interaction  energy  between  a  charge  distribution  U*  „ 

VA. 

(r.e.^,)  U  .  (r.8.<|>.)  and  another  charge  distribution  U* 

III  a  III  ||*jC»IU|| 

I  1  J  J  Zj 

( rje j ^un  i  m  (rj9 j# j )  where  it  is  understood  that  we  are  to  perform 

j  J 

a  spherical  average  of  the  charges  or  potentials.  To  do  the  integral. 


we  can  first  compute  the  quantity 
* 


Va/’jW  Va/W  © dVj 


(ID 


6 


HIJJI 
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This  is  the  electrostatic  potential  energy  (in  atomic  units)  of  the 
charge  distribution  of  the  jth  electron,  at  the  position  of  the  ith. 


Then,  we  multiply  this  by  the  charge  distribution  U*  „  (r.9.6.) 

ViV  1  1  1 


n.t.m  (r .9  -4> . )  of  the  ith  electron  and  integrate  over  dv . .  This  is 

llJC.lll  I 


exactly  what  Equation  4  says. 


In  addtion  to  this,  we  must  perform  a  spherical  average.  We  must 

take  the  total  charge  enclosed  in  a  spherical  shell  bounded  by  spheres 

of  radii  r.,  r  +  dr.,  and  replace  the  actual  charge  by  a  charge  dis- 
J  J  3 

tribution  in  which  this  same  charge  is  uniformly  distributed  through  the 

shell.  Now  the  charge  located  in  the  volume  bounded  by  r.,  r.  +  dr., 

y  J  J  J 

e.,  e.  +  de.,  <j>.,  <|> .  +  de .  is  the  volume  element  r.  sin  e.dr.d<|>.  times 


J  J  J  J  J  J  J  J  j  J 

the  product  of  the  wave  function  and  its  conjugate.  Hence,  the  charge 


located  in  the  spherical  shell  will  be  the  integral  of  this  quantity  over 
0  and  4  or 


I  (2i.+l)(£^-|mS,^j)!  2  /  x  f 

II  '■>'}’>*>  I 

rK  I 

L%  J  (c°seij  s,n6j  d8j 


ZlT  rH 

0  J  JO 


times  (-e)  if  we  are  using  ordinary  units.  The  integral  over  0.  is 

J 


given  by 


(Jl.+  |m-  |)l 
2  J  *j 


The  integral  over  $.  is  Zir .  When  we  combine  these  facts,  we  find 

J 

that  the  charge  located  in  the  shell  between  r.  and  r.  +  dr.  is 

J  J  J 


PniH.(rj^drj  t1mes  (‘e) 

J  J 


AFWAL-TR-83-1 1 30 


In  other  words,  as  far  as  the  spherical  average  is  concerned,  the  result 

2 

is  just  as  if  we  had  disregarded  the  functions  of  angles  and  used  Pn  { 
directly  as  the  radial  charge  density.  J  J 

To  find  the  electrostatic  potential  of  Equation  11,  we  now  wish  to 
find  the  potential  energy  of  an  electron  at  distance  r^  from  the  origin, 

in  the  presence  of  a  shell  of  radius  r.  and  of  the  charge  equal  to 

?  <1 

P  (r-)dr.,  using  atomic  units  in  which  the  potential  energy  of  inter- 

njfcj  J  J  2 

action  between  two  electrons  at  a  distance  r  is  By  elementary 

r 

principles  of  electrostatics  a  spherical  shell  of  charge  has  a  potential 
at  external  points  equal  to  that  of  a  point  charge  of  the  same  total 
charge  concentrated  at  the  origin,  and  at  points  inside  the  shell  the 
potential  is  constant,  equal  to  its  value  at  the  surface  of  the  shell. 

Then  the  potential  energy  of  the  electron  at  distance  r^  from  the  origin 

(in  the  presence  of  the  shell)  is 

h  pv1(rj>drj  ,f  r,- » rj  <’4> 

1  J  J 

h  Pn,*<rj>drj  1f  r1  «  rj  l'5) 

J  J  J 

The  total  potential  energy  of  the  electron  at  distance  r^  from  the  origin, 
in  the  presence  of  the  electron  whose  wave  function  is  given  by  Pn  (r^.), 
is  then  ^  ^ 


Potential  energy  =  ~  YQ 


where 


*  I?  (rj)dri  +  L. 
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The  notation  Yq^j.,  ~)  was  introduced  by  Hartree.  To  evaluate  the 
integral  (Equation  5), we  must  next  multiply  this  potential  energy  by  a 
spherically  averaged  charge  distribution  of  the  ith  electron  and  integrate 
over  the  coordinates  of  that  electron.  The  result,  which  we  shall  designate 
as  F°(n.n.;  n^),  is 


F-tV,;  Yj)  -  J  dr, 


p2  (r  \  p2 

Vi  1  Vj 


2 

rUT 


dr. dr. 

•  w 


07) 


where  r(B)  is  the  larger  of  the  two  variables  r . ,  r. ;  that  is  when 

r  >  r.,  r(B)  =  r  ,  and  when  r.  <  r.,  r(B)  *  r.. 

J  3  ^  J  J 

We  now  have  found  the  values  of  the  various  integrals  involved  in 
the  average  value  of  the  Hamiltonian,  and  Equation  6  can  be  rewritten 
in  the  form: 

(H)av  =  EfOKnflj)  +  E(pairs  i.j)  F°(ni^i;  n^.)  (18) 

3.  HARTREE  EQUATIONS  AS  DETERMINED  BY  THE  VARIATION  METHOD 

We  are  now  ready  to  carry  out  the  variation  of  the  one-electron 
orbitals  so  as  to  minimize  the  average  Hamiltonian  of  Equation  18. 

Since  we  wish  to  minimize  the  Hamiltonian  with  respect  to  variation  of 
each  one  of  the  radial  wave  functions,  we  can  handle  these  separately: 

(H)av  must  be  a  minimum  as  a  particular  function  P  is  varied.  The 

ni  i 

terms  of  Equation  18  which  depend  on  this  function  are 

(H^)av  =  Kn^)  +  (j  t  i)  F°(nij»i;  n^.)  (19) 


9 
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Let  us  vary  this  quantity  subject  to  the  condition  that  the  normalization 

integral  /*  P*  .  (r.)dr.  should  remain  equal  to  unity.  Using  Lagrangean 

0  "i^i  1  1 

multipliers,  A.  such  that  A *  =  -e„  .  where  e  .  is  the  eigenvalue  of 
i  i  n.£.  n.£.  3 

the  problem.  Then  1  11 


5<H1>4V  -  Vl,loP"i£1(ri)dr1  =0 


(20) 


then  computing  the  various  integrals  individually,  we  have: 


d2  MAi+l}  2z 
TT +  ~7~ 


dr 


=  2 


i  l 

2  A.  (A . .  i ) 


-  ij 


C  Vi  [• 


r.p  r r  t ‘i'Vi'  at  p  dr 

Jo  Vi  L'dSf  ^  Vi  dr>- 


(21) 


where  the  last  equation  follows  as  a  result  that  the  integral  is  real 
and  therefore  is  equal  to  its  own  conjugate.  Note  also  that 

C  Vi  <4> « v,*i  ■  &  «v)  V,  l  - 


Vi'^-i  1dri*(^5Vi)Pl"iti 

5Vi  l0  +  K  Vi  6W 


0  "  dri  ^Pn1Ai^ 


since  the  boundary  conditions  for  the  extremum  at  0  and  «  are  zero  and 

P  (r)  vanishes  at  0 
"i  i 


oo  so 


”Pn  l  (ar>  SPn  t  dr1  ■  f  A  <p„  1  >  SPn  l  dr1 

Jq  nini  ar.  n^  i  J0  drf  Vi  n.x,.  l 


10 


a;.  «jf>*  *f *‘v  . «'  ■'  •\iV‘  .y  ./.j 


and 
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where  the  minus  sign  was  removed  for  convenience. 


Next,  we  have 

SF°<"iV  njV  - 2 1  f/Vi  (ri)Vi(r<> 

pnJij(rj>  m  drtdrj 

’  2  j!  SPnil1<''l>  V,(ri)  ^V0  (njV  drf 


(22) 


Last  we  have 


6f  pn  i  dr<  s  2  I  6Pn  9  p„  9  dri 

J  o  Vi  '  J 0  Vi  Vi  1 


(23) 


Then  we  insert  these  values  into  Equation  2,  as 

2 

JQ  ni*i  .  L_  dr‘ 

n.Jt. 


d2  ^  W»-l  ^  2z 


f  6Pn  t  (r^  \-^  +  - V 
J0  Vi  i  L  dr^  ri 

+  K1  t  j)  ^Y0  (n.ly  -^) 

”  Pni*idr1  ’  ° 


(24) 


The  integral  (Equation  24)  must  be  zero  for  any  value  of  the 
variation  function  P  (r.)  which  means  that  the  remaining  part  of  the 

Vi  1 

integrand  must  be  zero.  Hence,  we  have 


‘  ^7  + 
L  drf 


d2  A  Wl*  2z 


V(Jil  VY°  'Vr^] 


(25) 
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Equation  25  is  the  fundamental  equation  of  the  Hartree  method. 

It  is  the  radial  wave  equation  for  an  electron  moving  in  a  spheri¬ 
cal  potential  produced  by  the  nuclear  change  of  z  units  and  by  the 
spherically  average  charge  distributions  of  all  other  electrons.  In 
other  words,  we  are  precisely  led  to  the  set  of  equations  which  would 
be  derived  from  Hartree's  original,  intuitively  derived  postulate  of  the 
self-consistent  field. 

Hartree's  equations  are  of  the  form  called  integrodifferential 
equations.  The  only  practical  method  of  solving  the  set  of  simultaneous 
equations  for  the  occupied  orbitals  is  a  method  of  iteration  or  successive 
approximations.  As  for  the  actual  method  used  in  obtaining  the  solution 
the  numerical  method  will  be  discussed  in  detail  in  a  later  section. 


AniAi.-Ta-33-U33 


1 


SECTION  IV 


FINITE-DIFFERENCE  APPROACH  TO  SCHRODINGER'S 
EQUATION  FOR  THE  SPHERICAL  PROBLEM 


In  the  case  oT  spherical  symmetry,  we  can  make  separation  of 
variables  to  Schrodinger's  equation  in  spherical  coordinates,  writing 
U.  as  the  product  of  a  function  R^  of  the  radial  vector  r  and  a  function 


of  Yim(0»<(>).  where  e  and  $  are  the  angles  in  spherical  coordinates  and 


i  and  m  are  quantum  numbers,  then 


(-!  a  -±; 


1  3  U 

+  -x - x - V 

r^sin  9  W 


Substituting  R^Y^te ,<(•)  for  in  Equation  26,  we  have 


^&PS  +  (‘rv)r2Ri]- 

fsjn9  _i  !!!»1 

r n9  88  J  ^9 


i  r  i  i_ 

{.  LsinQ  39 


Since  the  LHS  is  a  function  of  r  only  and  RHS  is  a  function  of  0 


and  4  only,  each  of  these  is  a  constant,  fc(t+l),  so 


1  X  fcln  0  !I*2L|  +  —L 

sinO  36  (s1n  9  “WJ  “1 


0  W 


+  =  0 


(■*  2M-.  - '  - 


A(A+1 


Ri  =  0 
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Si 


%'S 

i\ 


m 

$ 


ft 


V 

,"*‘4 

>- 


The  first  step  in  the  treatment  of  Equation  29  is  to  replace  R.  by 
P  1 

a  function  so  that  P^  =  rR^ ,  then  Equation  29  is  converted  as  follows: 

7‘2ra?1  +  ^i+l*i-v-M715Ri'0  (30) 

2  d2R  .  dR.  r-  -j  2 

r  — T  +  2r  dT  +  Lei  "  V]  r  Ri  '  A(£+1)Ri  =  0  (31 ) 


R  a  _1 

Ri  r 


dR.  dPi 

3r  rdr  7 
r 

d2Ri  d2P. 


dP..  2dP. 


r  dr  r  dr  r  dr  r  dr 


Substituting  Equations  33  and  34  into  Equation  31,  we  have: 


/■rd  P.  2dP.  2dP.>  rdP.  P.-t 

(l3L--a?Ltnlr!j+*r  ^- ?]♦«,*«  Pi 


dP.  P..' 


*U+1) 


r  Pi  "  °  (35) 


divided  by  r  we  have: 


r-  *  [e.  -V  - 


P1  =° 


Letting 


we  have 


(Ha)  g(r)=  -[e.  -  V  - 

L-  r 


s  9 (r)  P. (r) 
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The  most  practical  way  to  handle  Equation  37  is  numerical  inte¬ 
gration. 

1 .  ITERATIVE  SOLUTION  PROCEDURE 

Suppose  we  have  found  values  of  the  function  P  at  a  set  of  equal - 

spaced  values  of  r,  namely,  Pn  ^n-l  ’  pn  differing  by  increments  of  h 

of  the  variable  r.  Let  us  suppose  we  know  the  expansion  of  P  in  power 

series  around  the  value  r  of  r  involved  in  P  .  We  then  have 

n  n 


^  _ 1 »  .  W3  lA  s 


•Vtfi  -  **■  2  *5r**;‘  1 5r%*  *  5r  "i*  ^  ft  p:  +  —  (38) 

where  P^,  P^1'  etc.  are  successive  derivatives  of  P  with  respect  to  r. 


computed  at  r  .  We  then  note  that 


2h2  ^  2h4  Biv  .  2h6  nvi 


J  +  P  =  OD  +  *■"  D  4.  c-n  D,v  j.  c‘u  DVI 

n+l  *  n-l  2Pn  ?r  Pn  +  4T  Pn  ITT  Pn 


Pn+1  *  2P„  *  p„-l  “  ♦  W  '1*  * 


+  *  *  * 


(39) 


(40) 


since 


Pn*  =  9n^Pn^‘  From  E<luat*on  37 »  we  have 


4  , 


Vi  -  (2  *  g/)p„  -  p„.,  ♦  ^  C  + 


(41) 


If  the  intervals  of  h  are  small  enough  so  that  the  fourth  order 


h4p u 

n 


term  — ^ 0,  this  allows  us  to  compute  the  next  entry  Pn+1  in 


the 


table  if  P„  and  P^  .  are  known, 
n  n- 1 


Generally,  however,  this  approximation  is  not  accurate  enough.  It 
is  therefore  customary  to  use  a  modification  of  the  procedure  suggested 
by  Noumerov,  which  Is  more  accurate.  We  define  a  quantity  y^  by  the 
definition: 

9  I  l 

y  -  P"  Z3 

n  ^  (42) 


15 


AFWAL-TR-83-1 1 30 


since 


h2  .."  ,  h3  .  h4  ,.1v 


Vr».‘l«,»*TIn*T»»  +24yn+"' 


2  1 1 
JTP_  ' 


piii  p 

.  »  P„  l  h2 


h2pivs 


4  K"- ^  Kv- ¥] 


C  nVll  g  Dvm 

+  h  fpv  .  _n_l+  fpvi  _  In _ V  ... 

“  JI5  (Pn  12  J  27  [  n  12  J" 


i  o  •  •  i  ii  o  iii  r  i  -r 

Vl  *  p„  £  hpr.  +  h  pn  llT+?)jlh  P„  [rf  1  I] 


1,4  C  (il  *  ?s]  +  1,5  Pn  (t?  £  raj)  +  "6  Pn1  ((li)(24)+sr] 


+5  .2-''  h3  h5  „v  h6  Dv1  ... 


Vl  *  p„  *  hP„  it  h  P„  £f?pn  *  T55  n  ~  180  n 


from  which 


y  ,  +  y  i  =  2P  +  —  h2  P ' '  pvi 

yn+l  yn-l  cvn  +  6  "  Pn  WS  Pn 


or  for 


P  =  y  + 
n 


P  1 1 

h  P. 


we  have 


V,  -  2*n  £  «v,  ■  h2  p;  -^p*' 


i*  »•  »  *  .'p***! Rk • W 

O.  jl*. 
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u  r  * 

showing  that  the  fourth  term  has  disappeared.  Using  - 5-  =  g(r)  P.(r)„, 

h^P' 1  2  ^  2  1 

we  n,ay  rewrite  Vn  ’  =  P,  •  9„P„  "  Pe  ( 1  '  T2  !>„)• 

_  o  r1  ^ n  "  ^  2y -i  .6 

Vl  -  *»  *  Vl  *  h  [-^2 — ^  -  W  C  t«) 

yn+1  ■  2yn  •  Vl  +  [■  7 g%  -  12yn]  *  5W  Pn’  (50) 


y  =  [2  +  ^ - -]  y  -  y  .  pVi  , 

n+1  (.  pj  yr\  yn-l  240  n 


1  -  g  h 
yn 


am 


Here  the  sixth-order  term  is  small  enough  to  neglect  in  actual  cases, 
so  Equation  51  is  used  to  carry  out  the  Noumerov  integration.  First, 
this  is  done  by  constructing  a  table  of  values  of  the  quantities 
g  h2 

from  equation  g(r)  =  -  [ei  -  v  -  ^^-1.  Then,  we  carry  through 

*-  r  J 

two  columns,  one  of  y  ,  the  other  of  P  .  From  the  entries  in  the  y  table, 

n  n  n 

up  to  y  itself,  we  can  use  Equation  51  to  find  y  +1  and  from  equation 
n  9  n 

/  ~  i_£  \ 


yn  =  Pn  {'  '  hr  |  we  Set  Pr 


The  iterative  procedure  has  been  followed  in  the  succeeding  software 


program. 


Once  P  ,  P  .  and  P  .  has  been  calculated,  then  Simpson's  rule  can 
n  n- i  n+i 

be  employed  to  calculate  the  integral.  This  is  proven  as  follows: 


i  (r-r  )‘ 

Since  P  =  Pn  +  (r-r  )  P*  +  — Jl- 
n  '  n'  n  2 


pi i  +  . . 

n 


^  * 1  ^  itm  K VJJ,  JJ  J,  JflW.  W  ■;  »>«.l » l ».,  I.L  1.WHJ  ■,'.  1.,  i*  *:.ij.^  ;r  '.t'.t. 
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then 


r  «+h 


r  n  ■•  r'n  "  r  •  (r-r„)‘  • n 


dr 


Vh 


Vh 

,3  ,, 


?hJ  11  r  1 '  I 

=  pn2h+2h_pn  ♦••*2h[_Pn*VPn  +"J 


(53) 


since  h  Pn  =  Pn+^  -  2Pn  +  P  p  then  we  find 


r +h 


fn 

J  r 


r„-h  Mr  ■  2h  lPn  +  -S 


IV^r1- jp„  +  V+ 


(54) 


h 

3 


Pn-1  +  4Pn  +  Pn+1 


Equation  54  is  a  simple  form  of  Simpson's  rule;  instead 


rr"+h  fVh 

ivh  Jvh 


P^dr  which  can  be  employed  to  /“  p2dr  *  N  normalizing. 


we  have 


2 

/“’fP-l  dr  =  1  so  Simpson's  rule  can  be  used  to  find  definite 

0  L/nJ 

integrals  or  indefinite  integrals  up  to  a  given  entry  in  the  table  and 
its  procedure  is  followed  to  find  the  normalization  factors  for  the 
radial  wave  function. 


(Quadratic  integrands  do  not  follow  Simpson's  rule.  This  is  shown  as 
follows. 

rr  +h 


for 


r 

i  r_  ■ 


P2dr 


18 
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P2  =  [p„  ♦  (r-rn)pj  *  p''  *  •  •]LPn+  (r'rnK 


(r-rS  , 


)  .T| 

n__  pn 
n_ 


i  2(r-r  )  ii 

|l  x  H  n'1! 


■p„t2pn<r-rn>p:  +  — '.V  ‘"."i 


(r-r )3  ,  jj  (r-r  )4  ii2 

+  2 *3 —  pIp!1  +  TT—  Pi  + 


2  '  n  n 


integrating,  we  then  have 


P2dr  -  Pn22h  .  4  P«%  *4  P’2  ^  P"‘  * 

Jr  -h 
n 

■  f  P  +  h2pnipn  +  *£I  *  "• 


since 


A"  *  Pn+1  '  2Pn  *  Pn-1 


* 

3|l  (Ipn2  +  P  P  ,  -  2P2  +  P  iP  +  h2P*  + 
T  L.  n  n+1  n  n-1  n  nj 


_2hD2,DP  +p  p+  h2Pl2  +  ** 
3  Kn  n  n+1  n-1  n  n 


>  t  2  fp 

n  3  l__n 


+  Pn+1  +  Pr 


for  h  <.<1 


letting 


2P  h  b'  hn  r  “ 

n  _  D.  ->  <1  p  +  p  .  +  P  ,  *: 

3  3  T  l_n  n+1  n-lj. 


not  Simpson's  rule. 
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However,  when  point  by  point  is  squared  and  then  parabolically  integrated, 
Simpson's  rule  is  applicable  and  the  properties  of  the  recursion  formula 
are  preserved. 

With  this  procedure  as  well  as  simple  arithmetic  ones,  the  following 
quantities  are  calculated: 

radial  electron  density  £P2 

number  of  electrons 

actual  electron  density 

These  quantities  contribute  to  the  calculation  of  the  self-consistent 
potential  in  the  next  section. 

The  previous  discussion  touched  on  the  determination  of  polynomial 
coefficients,  but  the  picture  is  incomplete  since  the  total  power  series 
expressions  have  not  been  discussed.  Furthermore,  initial  values  for 
coefficient  have  not  been  determined  up  to  this  point.  Therefore,  we 
will  discuss  this  now. 

We  can  write  the  quantity  g  of  Equation  58 


2z 

where  V  =  has  been  substituted. 

r 

We  express  P  as  a  power  series 

P  =  E(n)pnrn  (58) 


20 
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We  can  then  express  the  terms  of  the  series  expansion  in  the  radial 
equation 

d*“P  •> 

~ J*  9(r)P(r)  -  Z(n)(n-l}p  r"'2  =  Z(n,m)  V  d  m  ,rn‘2 
dr  n  m  n-m-d 

where 

£(n)(n-1)pnrn'2 

dr 

and 

g(r)P(r)  *  E(n,m)  VmPn_m_2rn"2 

If  we  define 

V_2  =  1(4+1) 


(59) 

(60) 

(61) 

(62) 


V_1  -  -2z  (63^ 

Vm  *  0  for  m  <  -2  (64) 

m 


we  can  then  equate  the  coefficients  of  equal  powers  of  r  in  Equations  60 
and  61  and  obtain  the  recurrence  relations 


Z(n)n(n-l)pnr 


,n-2 


J  2  v-  pn-m-2r 


n-2 


n  m 


m 


or 

n(n-l )  -  I  VmPn_m_2  ■  Um)  (65) 

From  these  recurrence  relations,  we  can  find  values  of  all  Pn’s;  once  we 
have  started  the  recurrence,  we  can  get  an  explicit  expression  for  the 
power  series  expansion  of  the  function  P. 
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Writing  out  Equation  65  explicitly 

n(n-l)pn  *  !(£+l)pn  -  2z  pn_>,  +  VQpn_2  +  V1Pn_3  *" 

Combining  terms  for  p  : 

3  *n 

-  »(W)]  P„  *  -2zpn.,  ♦  V0pn.2  *  V0pn.3  ♦ 

We  can  expect  from  our  experience  with  the  hydrogenic  problem  that 

the  first  non-vanishing  term  in  the  expansion  of  P  will  have  n  =  £+1 . 

For  this  term,  the  coefficient  n(n-l)  -  £(£+1)  is  zero,  so  that  Equation 

67  is  identically  satisfied  with  an  arbitrary  value  of  P  i ,  combined 

wi  th  P  =0  for  n  <  £.+1 . 
n 

The  quantity  P  ^  is  the  arbitrary  constant  that  will  eventually  be 
determined  to  normalize  the  wave  function.  We  can  write  the  successive 
equations  of  Equation  67  in  the  following  form: 

(0)P£+1  -  0  =>  p£+1  is  arbitrary 

(2it+2)Pjl+2  =  -2Zp£+1  for  n  =  £+1 
2(2£+3)p£+3  *  -2z  PA+2  +  V0  pA+1 
3(2JL+4)p£+4  *  “2z  Pa+3  ♦  V0Pa+2  ♦  V,p1+1 
4(2£+5)Pjl+5  -  -2z  p£+4  +  V0Pa+3  +  V1Pjl+2  +  VgP^ 

Since  P£+^  is  arbitrary  for  £  =  0  (IS  electron)  we  let  Pl  =  1.0000  and 
iterate  onward. 
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SECTION  V 

SYNOPSIS  OF  SLATER'S  APPROACH  TO  THE  CALCULATION  OF 
THE  ELECTRONIC  STRUCTURE  OF  ATOMS  AND  MOLECULES 


The  highlights  of  the  self-consistent  approach  are  stated  here  to 
clarify  any  problems  the  reader  may  incur  while  pursuing  through  the 
general  theory.  This  will  help  one  in  not  losing  sight  of  the  main 
thrust  of  the  presentation  as  well  as  not  losing  continuity  of  thought 
through  the  labyrinthian  paths  of  the  theory  as  presented  previously 
and  in  the  following  sections. 


1 .  THE  MANY  ELECTRON  HAMILTONIAN 


The  many  electron  Hamiltonian  forms  an  atomic  or  molecular  system  in 
atomic  units: 


zZ. 


-  f  7i  '  fza  r,“  *  rfj  r*j  +  a  jib  2  r 


ZaZb 


ab 


(67) 


Here  the  index  i  refers  to  the  electrons;  the  index  a  refers  to  the 

nuclei;  ria  is  the  distance  between  the  ith  electron  and  the  nucleus  a, 

r. .  is  the  distance  between  the  electrons  i  and  j,  and  r  .  is  the  distance 
i  J  ab 

between  nuclei  a  and  b.  We  also  assume  that  the  nuclei  are  at  rest. 


When  the  Hamiltonian  is  allowed  to  operate  on  the  product  wave 
function  U^ (1 ) . . .U^(N) ,  multiplied  by  the  conjugate  of  the  wave  function 
and  integrated  over  all  values  of  3N-dimensional  space  of  the  3N  electronic 
coordinates,  one  has  the  total  energy  of  the  atom  or  molecule.  Additionally, 
one  also  has  to  take  into  account  the  antisymmetry  of  the  wave  function 
when  the  coordinates  (and  spins)  of  any  two  electrons  are  intercharged  - 
this  is  the  Pauli  exclusion  principle.  It  is  therefore  necessary  to 
write  down  the  many  electron  wave  functions  not  in  the  form  of  a  product 
U^ (1 ) . . . UN(1 )  but  in  the  form  of  the  Slater  determinant 


-Wry 


u! (1 )  U-j (2)  ...  Uj  (N) 
U2(l)  U2( 2)  ...  U2(N) 

UN(1>  UN<2>  •••  V") 


(68) 
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One  then  can  vary  the  orbitals  U,  to  minimize  the  total  energy 
This  total  energy  has  the  form  9y 

<Ehf>  -  -  Z  n,  /U*  v2  u(dv  -  I  n,  /  u*  VNU(  dv 
X  1  u*(l)  U*{2)  (;i-)  Uj(l )U^ (2)  dv,dv2  *  E  X  ^ 


a?*b  ab 


the  Jte,S  TT"  “*  be  eXPU1"ed  "  f0ll0"S-  Us1"5  fact  that 
total  charge  density  of  all  electrons  is  normally  .  x  u*u. ,  where  the 

sedation  goes  oyer  all  the  occupied  orbitals  and  .here  L  minus  sign 
n  .cates  that  the  electrons  have  a  negative  charge,  it  has  been  found 
■"ore  convenient  to  assign  occupation  numbers  n.  to  the  states  such  that 
un.ty  represents  an  occupied  state  and  zero  an  empty  one.  Then  the  total 
electronic  charge  density  I,  ,  .  .  *  n.u^  and  where  the  nuclear  potential 
2Za  ' 

is  Vn.e  — ,  the  electronic  potential  is  V  (1)  =  /„(2)  (J_) dv  a„d 
a  r12  2 

6  msimsj ^  indicates  that  one  includes  only  those  pairs  of  occupied 
orbitals  1 , j  for  which  the  spin  quantum  numbers  m  or  m  are  the  samea.e 
pairs  of  spin  orbitals  with  parallel  spins.  Furthermore^  correction  ” 
for  the  potential  acting  on  the  ith  electron  of  position  1  must  be  taken 
1nt°  3CC0unt  (Vxi  "  'V2>¥2>  (J~)  dv2).  This  arises  because  the 
electron  in  one  orbital  does  not  act  on  itself;  to  see  this,  let  us  write 

Ve(1)  +  Vxi(1>  =  /  P(2)  dv2  +  /  U*(2)  U.(2)  (X*  dv2 

’  /I*  Ji  "J  “J(2)  V2)]  <£>  *2  (70, 

then 

jy,  "j  uj(2>  uj(2)j 


24 


AFWAL- TR-83-1 1 30 


I 


i 


.tv 


*3 


KS 


m 


which  represents  the  total  charge  density  at  position  2  from  (N-l)  elec¬ 
trons.  The  ith  electron  is  being  acted  upon  by  the  (N-l)  other  electrons 


The  most  interesting  term  in  the  Hamiltonian  is 

1  r  —  —  .  / _  _  \  f  /  2 


i  ^  Vj  I  V1'  V2)  <4>  V’>  V2>  IV, 


x  dv. 


(71) 


Using  the  simple  case  of  two  electrons  whose  antisymmetric  wave  function 
is  =  U  (1  )U  (2)  -  U.(1)U.  (2)  then  the  operator  average  is: 

*  J  v  * 


If  i p*  Oipdv]  dv2  =  //(U*(l)  U*(2)  -  U*(l)  U*(2)} 


0  (U.(l)  Uj(2)  -  [}.(])  Ui ( 2 ) )dv^  dv. 


*  //|uj(l)  U*(2)  6  U.(l)  Uj (2)  +  U*(l)  U*(2)  0  Uj(l)  ^(2) 


-  U^(l)  U*(2)  0  Ujd)  U .  ( 2 )  -  U*(l )  U*(2)  0  U.(l)  Uj(2)  (72) 


x  dv^  dv2 


In  the  two  negative  terms  we  see  that  an  exchange  in  the  positions  of  the 


two  electrons  occurs  for  the  operator  0.  In  the  above  case  0  =  - ,  a 


12 


non-differential  operator.  These  exchange  terms,  as  represented  by 
Equation  72  in  the  general  case,  can  be  considered  a  correction  to  the 
direct  Coulomb  interaction  energy  due  to  the  Pauli  exclusion  principle 


and  also  contain  the  simple  correction  term  V  .  arising  because  an  electron 


does  not  act  on  itself. 
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to  each  of  the  spin  orbitals  U. ,  the  resulting  one-electron  Schrodinger 
equation  is 

-72Ui  -  (*M  *  ».)  «,  -  l  4("sf"sj>  /[UJ(2)  <4’  U,U) 


(73) 


x  Uj  =  eiUi 


Equations  69  and  72  constitute  the  self-consistent  Hartree-Fock 
method.  This  will  be  explained  in  more  detail  in  the  next  section. 


2.  THE  X-ALPHA  METHOD 

The  x-alpha  method  is  an  attempt  to  simplify  the  average  energy 
expression  by  approximating  the  exchange  terms  in  a  certain  way.  The 
approximation  is  based  on  the  theory  of  the  free-electron  gas  and  will 
be  discussed  in  the  next  section.  However,  for  comparison  with  the 
Hartree-Fock  energy  (Equation  69)  the  total  energy  of  the  molecule  in 
the  xo  method  is  written: 

<Exa>  =  *E(i)ni  /  Ui  y2ui  dv 

*  2Za  i  ? 

-  Z(i)  n.  /  U .  E(a)  (— *-)  U.dv  +  l  p(1)p(2)  (/-)  dv,  dv2 

ii  rla  1  c  r12  c 

-  7a(^1/3  |{[-  p+0)^4/3  +  C-pHD]4/3}  dv, 

2Z  Zn  (74) 

+  £(pairs  aB,  a^B)  —— — 

raB 

Here  n,.  and  have  the  same  significance  as  discussed  earlier  as 
p,  Z  .  r,  ,  r17,  and  r  R.  The  quantities  [-  p+(l)]4/3  and  [-  pi(l)]4/3 
arise  from  the  product  of  the  quantity  [-  p+(l)]  '  or  [-  p  +  (l)]  and 
the  corresponding  charge  density.  The  value  of  the  exchange  term,  the 
term  proportional  to  a  in  Equation  74  is  based  on  the  theory  of  the 
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free-electron  gas.  The  arrows  (f)  and  (+)  denote  the  spin-up  and  spin- 
down  charge  densities,  respectively. 

From  Equations  69  and  74  one  can  see  that  the  xa  method  differs  from 
Hartree-Fock  only  in  the  exchange  term. 

3.  THE  EVALUATION  OF  THE  SELF-CONSISTENT  POTENTIAL 
FOR  AN  ATOM  IN  THE  Xa  METHOD 

Looking  at  the  various  potentials,  the  calculation  of 

V’>  ■  77  (75) 

is  obvious.  The  evaluation  of 


However,  for  a  spherical  distribution  p(r2)  of  charge,  we  can  write 
V  (1)  '  -T-  f  1  47ir l  p(r.)  dr.  +  P*  ~  4itr?  p(r,)  dr,' 

e  rl  JO  2  2  2  Jr,  r2  2  2  2  (77) 

This  expresses  the  familiar  fact  that  the  charge  within  a  sphere 

rl  2 

of  radius  r1 ,  equal  to  /Q  4*^  pfrg)  d^  exerts  a  potential  at  radius 

! 

r^  as  if  it  were  all  concentrated  at  the  center  of  a  sphere,  whereas  the 

charge  in  a  spherical  shell  of  radius  from  r2  to  r2  +  dr2  produces  a 

constant  potential  at  interior  points  equal  to  the  potential  just  outside 

the  shell.  The  integral  in  the  first  term  of  Equation  68  is  given  as 

rl  2 
-  /g  EP  dr. 

The  second  integral  can  be  immediately  computed  from  Simpson's 
rule,  after  constructing  a  table  of  functions. 
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The  next  potential  to  be  considered  is 


vx  ■  6“(Jr),/3  (-  p(D)'/3 


2  2 

The  charge  density  -p  =  zP  /4ir r  can  be  tabulated  so  that  V  .  is 

A  I 

essentially  computed. 


The  table  can  then  be  set  up  as 

r>  -Ve  Vx  VN  Vve  v  =  +  V.  +  V 


N  'e  'x 


This  table  will  give  values  of  total  potential  V  acting  on  the  electron 
and  its  Coulomb  part  (VN+Vg). 

4.  THE  TOTAL  ENERGY  OF  AN  ATOM 

This  involves  the  calculation  of  various  integrals  and  the  methods 
described  in  the  preceding  sections  allow  us  to  compute  them. 

The  one-electron  Schrodinger  equation  for  the  atomic  orbital  Ui ,  is 


where 


v2u(  --(«,+  V)U( 


v  -  VN  +  ve  ♦  vx 

v-a 

r 


ve  ‘  /p<r2><^>  dv2 

V.  ■  fa  (5r)1/3  [-  p(l)],/3 


From  the  previous  sections,  we  have  seen  how  to  compute  Vg  and  Vx  as 
functions  of  r;  the  computation  of  is  somewhat  trivial.  If  we  rewrite 
the  total  energy  <Exa>  in  terms  of  integrals  over  these  potentials. 
Equation  74  takes  the  form 
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<Exa<  =  -  E(t)ni  Jo*  V2u.dv  -  Eton,  J U*  VN  U, 


/U,  Ve  U.  dv-foEdln,  /U,VxU(dv  {g4) 

If  we  multiply  Equation  79  by  -U.  and  integrate,  we  have 

-  /IL  V2U.  dv  =  e.  +  /IL  V  U.  dv  +  /U*  V  U.  dv  +  /Ui  V  U.  dv 

11  1  1  n  1  161  (85] 

Thus,  we  see  that  if  the  U^'s  satisfy  Equation  67  we  can  reduce  the 

calculation  of  <E  >  to  the  evaluation  of  the  integrals  / U?  V„  U.  dv, 

*  i  N  i 

/U..  Vg  U.  dv  and  /Ui  Vx  U.  dv,  for  the  occupied  orbitals.  For  example, 
the  only  occupied  orbitals  tor  carbon  are  Is,  2s  and  2p,  therefore  this 
means  only  nine  integrals. 


Each  of  the  integrands  is  the  product  of  a  charge  density  U.  Ih  and 
a  function  VN,  Vg  or  Vx  of  r  only.  The  integrations  over  angle,  as  in 
the  problem  of  normalization,  lead  to  the  result  that  the  integrals  in 

r\  eo  P 

question  reduce  to  /g  Vn  dr,  /  Vg  dr,  and  so  forth.  We  can 

2  ® 

multiply  at  once  the  functions  P1$  by  VN  or  Vg  or  Vx,  and  carry  out  the 
integrations  over  r  by  Simpson's  Rule.  This  procedure  is  followed  in  the 
computer  software. 
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SECTION  VI 

THE  HARTREE-FOCK  METHOD  AND  ITS  SIMPLIFICATION 

The  Hartree  equations  were  obtained  by  varying  the  one-electron 
wave  functions  U.j(x),  ^(x)  ...  Un(x)  in  such  a  way  as  to  make  the 
energy  U*(x^)  ...  Un(xn)  H  U^(x^)  ...  Un(xn)  dxi  ...  dxn  an  extreme. 
This  was  derived  earlier. 


Similarly,  the  Hartree-Fock  equations,  which  take  into  account  the 
antisymmetry  properties  of  the  electron  wave  functions,  are  obtained  by 
varying  the  IKs  so  as  to  make  the  energy: 

1  r  ui (xi )  •••  u](xn)  ui^xi^  Mxn) 

—  i  — j j -  H  -  dx,  . . .  dx„ 

1  W  •••  u„<‘n>  W  •••  m*„)  1 

an  extremum,  where  in  this  latter  expression  the  U's  are  assumed  to  be 
functions  depending  on  coordinates  and  spins,  and  where  the  integrations 
over  the  dx's  are  interpreted  to  include  summing  over  the  spins.  The 
Hartree-Fock  equations  can  then  be  written  in  the  form: 


wv  ♦  [j,  K<*2>  w  w 

-  j,  [K<*2>  Ui<*2>  [4^=77] dx2  W  ■  W*1 


Here  H^  is  the  kinetic  energy  operator  for  the  electron  of  coordinate 
x.,  plus  its  potential  energy  in  the  field  of  all  nuclei:  * — - -  is 

1  4ireorl  2 

the  Coulomb  potential  energy  of  interaction  between  electrons  1  and  2 
expressed  in  MKS  units.  The  u^ 's,  as  before,  are  assumed  to  depend  on 
spin  as  well  as  coordinates,  and  the  integrations  over  dx2  include 
summation  over  spin,  so  that  the  exchange  terms,  the  last  ones  on  the 
LHS  of  Equation  86  automatically  vanish  unless  the  functions  U.  and 
correspond  to  spins  in  the  same  direction. 
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1.  INTERPRETATION  OF  THE  HARTREE-FOCK  EQUATIONS 

We  can  subdivide  the  total  density  of  all  electrons  into  two  parts 

p+  from  those  with  plus  spins,  p-  from  those  with  minus  spins;  the  two 

* 

together  add  to  the  quantity  -  ei(k  =  l,...,n)  U^(x)U^(x).  We  can  show 
that  the  Hartree-Fock  equation  (Equation  86)  for  a  wave  function  U.. 
which  happens  to  correspond  to  an  electron  with  a  plus  spin  is  an 
ordinary  Schrodinger  equation  for  an  electron  moving  in  a  perfectly 
conventional  potential  field.  This  field  is  calculated  by  electro¬ 
statics  from  all  the  nuclei,  and  from  a  distribution  of  electronic 
charge  consisting  of  the  whole  of  p -,  but  of  p+  corrected  by  removing 
from  the  immediate  vicinity  of  the  electron,  whose  wave  function  we  are 
finding,  a  correction  or  exchange  charge  density  whose  total  amount  is 
just  enough  to  equal  a  single  electronic  charge.  That  is,  this  corrected 
charge  distribution  equals  the  charge  of  n-1  electrons  as  it  should. 

The  exchange  charge  density  equals  just  p+  at  the  position  of  the  electron 
in  question,  gradually  falling  off  as  we  go  away  from  that  point.  We  can 
get  a  rough  idea  of  the  distance  in  which  it  has  fallen  to  a  small  value 
by  replacing  it  by  a  constant  density  p+  inside  a  sphere  of  radius  rQ, 
zero  outside  the  sphere.  We  have 

j  it  rQ  1  p+ 1  *  e  or 


The  situation  is  then  much  as  if  the  corrected  charge  density 
equaled  the  actual  total  electronic  charge  density  outside  the  sphere, 
with  p-  within  the  sphere;  this  is  the  Fermi  or  exchange  hole  surrounding 
the  electron  in  question,  consisting  of  a  deficiency  of  charge  of  the 
same  spin  as  the  electron  in  question. 

This  exchange  hole  is  clearly  different  for  wave  functions  of  the 
two  spins,  provided  p+  and  p-  are  unalike,  furthermore,  it  is  different 
for  each  dissimilar  wave  function  .  It  is  this  difference  which  leads 
to  the  complicated  form  of  the  Hartree-Fock  equation  (this  will  lead  us 
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to  the  simplification  of  using  a  sort  of  averaged  exchange  hole  for  all 
the  electrons).  The  exchange  holes  for  dissimilar  U.'s  of  the  same  spin 
will  only  show  small  differences.  We  have  already  seen  that  the  radius 
r^  which  we  obtained  by  assuming  a  hole  of  constant  density  depends  only 
on  p+  (for  a  plus  spin),  and  hence,  is  the  same  for  all  ’ s  of  that  spin. 

It  is  this  small  dependence  on  which  will  make  it  reasonable  to 
simplify  the  Hartree-Fock  equations  by  using  an  averaged  exchange  charge. 

To  agree  with  the  above  description,  we  expect  that  the  exchange 

charge  density  at  point  x2,  producing  a  field  acting  on  the  electron  ait 

,  whose  wave  function  U^(x^)  we  are  determining  by  the  H-F  equation 

(Equation  86),  will  have  a  volume  integral  equation  to  -p  and  will  be 

n  * 

equal  when  x?  approaches  x,  to  the  quantity  -e  l  U.  (x, )  U.(x,). 

k=]  k  i  k  1 

To  show  this  we  rewrite  Equation  86  in  the  equivalent  form 


HjU.^)  + 


J,  K  WV  (; 


4”Vl2 


dx2  Ui(x]) 


jJl  ^MX1*W  WW  4Tve0r12)dx2]  W 


Ui(x])  Ui(x1) 


=  WX1> 


(88) 


The  exchange  term  now  appears  as  a  product  of  a  function  U. (x.j) 
times  the  function  U^x^).  This  exchange  potential  energy  is  the 
potential  energy,  at  the  position  of  the  first  electron,  of  the  exchange 
charge  density 

_e  l 

k*1  Vxl>  (89) 

located  at  the  position  x2  of  the  second  electron. 
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We  note  that  the  exchange  charge  density  depends  on  the  position  of 

the  first  electron,  as  well  as  the  second,  and  also  on  the  quantum  state 

i  in  which  the  first  electron  is  located.  Furthermore,  the  total  charge 

is  that  of  a  single  electron;  also,  as  approaches  x^ ,  we  see  at  once 

that  the  exchange  charge  density  approaches  the  correct  value 
n 

-e  T.  uJ{x1)UR(x1).  Thus,  we  see  that  the  exchange  term  (Equation  89) 

satisfies  all  the  conditions  necessary  to  justify  our  qualitative  discussion 
of  its  behavior. 


AVERAGED  EXCHANGED  CHARGE 


The  probability  that  an  electron  x-j  should  be  in  the  i  state  is 
evident 


Uj (X] )Ui (xi ) 


(90) 


We  can  use  this  quantity  as  a  weighting  factor  to  weigh  the  exchange 
charge  density  (Reference  3).  When  we  do  this,  we  find  as  the  average 
exchange  density,  the  quantity 


n  n 


j WWWUj 


(x9) 


jJ,  Uj(xl)UJ(xl) 


(91) 


Using  this  average  exchange  charge  density,  we  obtain  in  place  of  the 


Hartree-Fock  equations  the  following  Schrodinger  equations  for  the  U^'s: 


H1Ui(x1)  +  ■ 


J1  ^Vx2)lVx2* 


4lTe0rl  2 


dx. 


n  n 


j k^  ^Uj(xl)Uk(x2)Uk(xl)Uj(x2) 


4lT€0r12 


dx 


=  EiUi(x1) 


W 


(92) 
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The  wave  functions  and  energy  values  ,  as  determined  from  these 
equations,  will  not  be  quite  so  accurate  as  those  determined  from  the 
H-F  equations,  but  they  will  at  least  be  much  better  than  those  found 
from  the  Hartree  equations. 

One  aspect  of  Equation  92  that  will  prove  to  be  important  is  that 
the  functions  U.  are  all  orthogonal  to  each  other. 

3.  EXCHANGE  CHARGE  FOR  THE  FREE-ELECTRON  CASE 

The  calculations  of  exchange  charge  and  exchange  potential  which  we 
have  been  describing  can  be  carried  out  exactly  for  the  case  of  the  free- 
electron  gas.  In  turn,  the  results  for  the  free-electron  gas  can  be 
directly  applied  to  approximate  the  exchange  potential,  with  an  expression 
which  is  much  simpler  than  that  appearing  in  Equation  92. 

With  this  in  mind,  we  now  develop  the  theory  of  the  free-electron 


Let  us  consider  a  free-electron  gas  with  N  electrons  in  a  cube  of 
side  L.  The  Fermi-Dirac  distribution  function  is 

f  n  (e-6r)kt  1 

fF'D(e)  =  [e  F  +  1]  1  g(e)  (93) 

where  fp  is  a  parameter  characteristic  of  the  electron  gas  and  g(e)  is 
the  total  density  of  states  per  unit  volume.  To  calculate  the  latter 
we  must  confine  ourselves  to  the  points  k  which  are  contained  within 
the  surface  defined  by 

kx2  +  ky2  +  kz2-^=0  (g4) 

We  now  wish  to  determine  the  number  of  allowed  energy  levels  in  terms  of 

the  number  of  different  sets  of  positive  integers  (m  ,  m  ,  m  )  which  give 

x  y  z 

rise  to  independent  solutions  of  the  wave  equation  for  free  electrons  in 
a  box: 


V  v  +  (E-V)t^  *  0,  where 


V  =  0  (0  <  (x,y,z)  <  L)  and 

V  =  oo  elsewhere. 
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Since  the  solution  of  Equation  95  is  ip 

sin  kxx  sink^y  sink2z  where 


W"z 


(x.y.z)  =  (|)1/2  x 


fk 

n 

X 

X 

k 

n 

y 

= 

y 

k 

n 

zJ 

zJ 

j-,  the  energy  eigenvalue  which 


2  2 

corresponds  to  the  state  (n  n  n)  has  the  form:  E(k)  =  =  e 

x  y  z 


Tjl  (1)  If-2  *  ny2  +  ,  V 

5m  V 


Since  n  *  -£-,  then  Equation  94  is 

X  TT 


n!  +  n!  +  n!‘T76s0-  (V  ls  the  volume  L3) 

X  “  Z  TtC 


(96) 


Thus,  the  number  of  states  may  be  taken  to  be  the  volume  of  the  first 
octant  of  the  sphere  in  n  space  which  is  bounded  according  to  the  surface 
(Equation  96). 

If  we  denote  the  total  number  of  states  at  the  energy  e  by  v(e), 
then  using  Equation  96  we  have 

v(e)  *1  %n3- 


[2m  L2  ’1 3/2 

1 

r2m) 

1?^6J 

= 

6tt 

ffl 

% 

1 


V  e 


3/2 


p(e)aM|l^[^3/2Vel/2 
d  e  4ir 


(97) 


(98) 


2  2  2  2 

where  n  =  nx  +  ny  +  nz  and  p(e)  is  the  density  of  states  of  a  given 
spin.  Then  the  total  density  of  states  per  unit  volume  is 

1 3/2 


9(e)  =  i  fe]  "  e1/2 

2/  kZJ 


(99) 
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and  the  average  value  of  any  physical  property  Q(r)  associated  with  the 
one-electron  energy  e  (per  unit  volume)  is 


Q(ej  e1/2de 
exp(e-er/KT] 


(100) 


The  concentration  of  free  electrons  according  to  Equation  93  is: 

i  r->m i3/2  r“  e^2de 


1  [2m }6/d  r 

oJ-  J  n 


0  exp(  F)  +  1 
KT 


(101) 


at  T  =  0°K, 


expi-j^)  +  1 


(102) 


therefore,  one  is  for  e  <  £p  and  zero  for  e  >  Ep  so  that  the  concentration 
for  a  highly  degenerate  electron  gas  (i.e.,  one  for  which  kT  «  Ep)  is 


Thus, 


n  3 

1 

2m] 3/2 

f6F  J/2 

1 

i7 

y 

e  de  = 

J0 

:  T 
3n 

eF 

.7 

'  2i 

(3n2n)2^2 

II 

rJ=£> 

II 

M— 1 

II 

CVJ  U_ 
> 

B 

-or 

PF 

h2 

[3n2n]2/3 

1  h2  fj 

InV 2/3 

2m 

*  S 

l  .  J 

,7 ' 2m  l« 

hr 

(103) 


(104) 


Therefore,  the  electrons  (for  the  degenerate  case)  occupy  energy 
values  with  uniform  density  in  momentum  space,  out  to  a  level  whose 

PF  ,,J/3 

energy  is  ^  ,  corresponding  to  a  maximum  momentum  Pp  =  h 


The  DeBroglie  wavelength 
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associated  with  this  maximum  momentum  is  related  to  the  radius  rQ  of  the 
exchange  hole.  We  notice  that  for  the  free-electron  gas,  the  |p+|  which 
appeared  there  equals  so  that 


(106) 


Finally,  in  order  to  calculate  the  exchange  term,  we  can,  with 
negligible  loss  of  accuracy,  describe  each  state  by  the  plane  wave 

4»K(r)  =  -1  e ^  where  K  =  ( V ny,nz)  and  Vny  and  nz 

/?" 


(107) 


are  integers.  The  exchange  integral,  defined  from  the  exchange  term  in 
the  H-F  Hamiltonian  (Equation  88): 

(Xj)  /W°k*x2*  (4ire0r12)  Uk*xl  )Ui  ^x2)dx2 


6ik 


UjUj )Ui(x1 ) 


e(K.K')  = 


4Tre0|rrr;|*k(f2^kl(fl)dr2 

^k(rl)^k(rl) 


j  irt-k1')  o  (r,-r2)  g2 


~7  \  e 


*reoirrr2 


dr-jdr2 


(108) 


4  e 


lK-K'1  2 


The  minus  sign  indicates  a  decrease  in  the  electron-electron  interaction 
between  electrons  of  the  same  spin  because  of  the  Pauli  principle. 
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If  K  is  the  wave  vector  for  a  given  electron,  then  K  is  an  index 
denoting  wave  vectors  for  all  other  electrons  in  the  system.  From 
Equation  88,  the  total  exchange  interaction  acting  on  a  single  electron 
is  then  obtained  by  summing  over  states  K's.  The  summation  can  be 
replaced  by  an  integration,  since 


/<*S  * 


T  1  up 

(2*h)3  (2tt) 


/  dK 


(109) 


Using  Equations  109  to  108,  the  total  exchange  energy  acting  on  a 
single  electron  is  thus 


£  e(K,K')  =  -—-2  fK|:  —I"? 

K'  47re^2i)T  > 0  |t-K*r 


2  r2ir  fir  r  Kc 

— a  d0  sinOdQ  h 

(4-nefiiif  J0  J0  J0 


K,2dK' 


K2+K'2  -  2KK'cos0 


*  £(^)  (110) 


where  Kp  is  the  Fermi  wave  vector,  i.e.. 


mer 


(111) 


Evaluation  of  the  integral  is  straightforward  but  lengthy  and 


gives  e(K) 


«%  f  K2  -  K2  |Kf-K|  , 

T2W7^Tre0  KKp  Kp+K  J 


(112) 


The  exchange  potential  energy  can  be  conveniently  stated  in  terms 

of  the  ratio  n  *  •§-  of  the  magnitude  of  the  momentum  of  the  electron 
F 

to  the  maximum  momentum  corresponding  to  the  top  of  the  Fermi  distri¬ 
bution. 
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The  exchange  potential  energy  =  - 


F(n ) 


where 


(113) 


(114) 


The  function  F(n )  is  shown  in  Figure  1.  It  goes  from  unity  when  n  =  0 
for  an  electron  of  zero  energy  to  j  when  n  =  1 .  at  the  top  of  the  Fermi 
distri bution. 


We  see  that  this  exchange  potential  energy  is  of  the  form  which  we 
expect.  If  we  had  a  sphere  of  radius  rQ,  filled  with  uniform  charge 
density  |p+|  *  the  potential  energy  of  an  electronic  charge  at  the 

center  of  the  sphere  would  be  -  |  while  the  value  of  Equation 

113  is  -1.54  {  r  c-t—  )  at  the  bottom  of  the  Fermi  band  and  half  this  value 


-1 54(^k) 


at  the  top.  Thus,  this  simple  model  of  an  exchange  hole  of  constant 
charge  density  gives  not  only  a  qualitatively  correct  value  for  the  exchange 
potential  but  also  a  rather  accurate  quantitative  value;  and  the  extreme 
difference  between  top  and  bottom  of  this  band  corresponds  only  to  a 
factor  of  2  in  the  exchange  potential.  If  we  average  overall  wave 
functions,  we  find  that  the  properly  weighted  average  of  F(n)  is  3/4. 


Thus  the  exchange  potential  energy  of  the  averaged  exchange  charge 


KU!U.kU.l  Hi  11U1BLU  W  A  "U'i1  A'AV. 
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4.  THE  USE  OF  THE  FREE-ELECTRON  APPROXIMATION  FOR 

THE  EXCHANGE  POTENTIAL 

It  is  dear  that  the  exchange  charge  density  (Equation  91)  and  the 
corresponding  potential  appearing  in  Equation  92  must  depend  on  the 
density  of  electronic  charge,  but  not  greatly  on  anything  else.  Thus, 
in  no  case  will  we  expect  it  to  be  very  different  than  what  we  should 
have  in  a  free-electron  gas  of  the  same  charge  density.  For  a  non-spin 
polarized  system,  in  which  p+  =  <r  *  we  may  approximate  the  averaged 
exchange  potential  by  what  we  should  have  in  a  free-electron  gas  of  the 
same  density.  Thus,  combining  Equations  87  and  115,  we  have 


eav(x)  = 


3  l6\ 
4-  (7) 


3/2 


(116) 


where  p(x)  is  the  local  electron  number  density  at  point  x.  If  we  recall 
that  this  is  z  Uk(x)Uk(x)  we  finally  have  our  simplified  Schrodinger 

equation  for  the  one-electron  function  ,  to  replace  Equation  92 


H^.U)  +  I  /Uk(x2)Uk(x2) 
K 


4ire0r 


d 


dx. 


(117) 


-3  [I-  jj  uk<xi>vv],/3  <M*1>-Wxi) 

This  equation  is  in  practice  a  simple  one  to  apply.  We  solve  it 
for  each  of  the  wave  functions  U.. ,  then  find  the  total  charge  density 
arising  from  all  these  wave  functions  and  then  calculate  the  potential 
energy,  including  the  exchange  term,  to  go  into  Equation  117,  so  as  to 
check  the  self-consistency  of  the  solution.  We  can  change  to  atomic 

units  by  setting  7^ —  to  2. 

4,re0 

This  method  will  facilitate  an  explicit  self-consistent  formulation 
for  the  case  of  atoms  as  follows.  Let  the  electrostatic  potential  of  the 
nucleus,  and  of  all  electrons,  at  distance  r  from  the  nucleus  of  the 


41 


'0 


m 


,  is 


y 
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7  (p\g 

spherical  atom,  be  _£___.  Then  the  charge  density  is  given  by  Poisson's 


equation  as 


4ireor 


c  =  -e0 9  <?%> 


018) 


When  we  express  the  Laplacian  in  spherical  coordinates,  this  gives 


o  ■  -  *>  <f>  Tf 

dr 


(119) 


The  exchange  potential  energy  becomes 


- -3  ia-  F^fkir/3 

* e0  32/  r  dr  * 


(120) 


and  finally,  the  total  potential  energy,  for  use  in  the  Schrodinger 
equation  for  (x . )  is 


d_  r  +  3  (-a,,’/*  p  %)m\ 

"0  L  P  32ir  (  dr  >  J 


(121) 


To  carry  out  a  self-consistent  solution  for  an  atom  using  this 


method,  we  find  Zp(r)  such  that  the  wave  functions  U. ,  determined  from 


a  single  Schrodinger  equation  using  the  potential  energy  (Equation  121) 


added  to  give  a  charge  density  which  would  lead  by  Poisson's  equation 


back  to  a  potential  energy  (— — £). 

AttEq 


When  this  procedure  is  carried  out  for  the  atoms  of  the  various 
elements,  it  is  found  that  the  ground-state  energy  of  an  atom  is  con¬ 
sistently  greater  than  that  obtained  in  the  Hartree-Fock  approximation. 
Since  the  H-F  equations  satisfy  a  variational  principle,  they  must  give 
an  upper  bound  to  the  energy  of  the  ground  state.  To  improve  the 
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agreement  with  <EU  ■>,  the  above  approximate  treatment  of  exchange  is 
modified  by  Including  a  parameter  n  to  represent  the  effects  of  the 
non- free  electron  nature  of  the  system.  The  value  of  a  is  determined 
by  setting  the  atomic  ground-state  energy  equal  to  the  Hartree-Fock 
value.  This  results  in  a  values  ranging  from  about  0.75  for  the  lighter 
elements  to  about  0.69  for  the  heavier  elements.* 

If  one  also  takes  into  account  spin  polarization  of  the  electrons, 
then  the  final  form  of  Slater's  x^  local  approximation  to  the  exchange 
potential  at  point  x  for  an  electron  of  spin  up  or  down  is,  in  atomic 
units, 

V*> = 

(122) 

where  p(x)  denotes  the  appropriate  spin-up  or  spin-down  electron  number 
density  at  x. 

Equation  92  for  the  one-electron  orbitals  U.  and  energies  E^  of 
an  atom  or  molecule  becomes 

[H-j  +  Vg(Xj )  +  Vxa(x1)]  U^(x-j)  s  E.j  (x.j )  ( 1 23 ) 

with  Ve  given  by  Equation  96  and  VXq  by  Equation  122,  and  with  total 
energy  <EXa>  as  given  in  Equation  94. 


*The  a  coefficient  can  also  be  chosen  so  as  to  satisfy  the  virial  theorem, 
<T>  =  -  1/2  <V>,  where,  <T>  and  <V>  are  respectively  the  average  kinetic 
and  potential  energies.  This  theorem  must  hold  for  the  true  ground  state 
of  any  system,  such  as  an  atom,  in  which  the  interactions  are  all  Coulombi 
The  values  of  a  determined  in  this  way  are  essentially  the  same  as  those 
determined  by  the  first  method. 
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SECTION  VII 

SIMPLE  SCATTERING  THEORY 

We  shall  consider  the  scattering  of  a  stream  of  charge  particles  by 
a  small  spherically  symmetrical  region  in  which  the  potential  energy  is 
different  from  zero. 


In  experiments  on  the  scattering  of  a  beam  of  particles,  one  measures 
the  number  of  scattered  particles  falling  per  unit  time  on  an  area  ds 
placed  at  a  distance  r  from  the  scattering  atoms.  The  number  of  particles 
falling  on  ds  is  proportional  to  the  area  ds  and  inversely  to  the  square 
of  the  distance  r.  We  shall  refer  to  the  particles  which  hit  ds  as 
"scattered  through  an  angle  e  into  the  solid  angle  dw."  This  number  of 
particles  scattered  into  dw  is  also  proportional  to  the  current  per  unit 
area  in  the  incident  beam.  Suppose  that  N  particles  cross  unit  area  per 
unit  time  in  the  incident  beam.  Let  the  number  of  particles  scattered 
per  unit  time  through  an  angle  0  into  the  solid  angle  dw  be: 

NI(Q)du  (124) 

Then  1(e)  is  the  quantity  we  wish  to  calculate.  1(e)  dw  has  the  dimensions 
of  area  and  is  called  the  effective  cross-section  for  scattering  into  the 
solid  angle  dw  or  as  the  differential  cross-section. 


Let  (x,y,z)  denote  Cartesian  coordinates  of  the  electron  and  (r,e,<|>) 
its  spherical  polar  coordinates.  We  shall  suppose  the  atom  to  be  situated 
at  the  origin,  and  the  potential  energy  of  an  electron  distance  r  from 
the  origin  to  be  V(r).  We  shall  also  assume  that  V(r)  tends  to  zero 
faster  than  1/r;  the  case  of  Coulomb  scattering  will  be  addressed  later. 

We  represent  the  stream  of  electrons  by  the  plane  wave  exp(iKz)  where  K 
is  equal  to  v-  =  This  wave  represents  a  density  of  electrons  per 

unit  volume,  and  therefore  a  flow  of  v  electrons  across  unit  area  per  unit 
time.  The  wave  will  be  scattered  by  the  atom,  the  amplitude  of  the 
scattered  wave  at  the  point  (r,  e,  $)  being 


f(Q)4iKr 

r 


(125) 
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that  is,  the  plane  wave  will  be  scattered  as  a  "wavelet"  by  Huygen's 

1  Kr 

principles  or  - -  but  whose  direction  is  modified  through  a  given  angle 

0.  Our  problem  is  to  find  f (0 ) .  From  it  we  can  deduce  the  number 
scattered  into  a  given  solid  angle  per  unit  time.  The  number  of  electrons 
in  the  scattered  wave  crossing  an  element  of  area  ds  at  the  point 
(r,  0,  <j>)  is  £  ds  |f(e)|2  per  unit  time;  and  therefore,  if  the  incident 
beam  is  such  that  one  electron  falls  on  unit  area  per  unit  time,  the 

number  I(e)dui  scattered  into  a  given  solid  angle  dui  per  unit  time  is 

2  7 

equal  to  | f ( 0 ) I  dm.  Therefore,  1(0)  =  |f(e)|  .  The  number  of  particles 

2 

scattered  between  angles  0  and  de  +  0  is  |  f (0 )  |  2i;sin6de,  Our  problem 
is  to  find  a  solution  ¥  of  the  wave  equation  which,  at  a  large  distance 
from  the  atom,  represents  an  incident  wave  and  a  scattered  wave,  or  for 
large  r. 


¥  ~  elKz  + 


(126) 


The  wave  equation  satisfied  by  v  may  be  written  as 


V2V  +  [K2  -  U(r)]¥  =  0 


where 


K  =  Sf;  U(r)  =  JrV(r) 


(127) 


Before  considering  the  solution  of  Equation  127,  we  require  a 

i  Kz 

certain  expansion  in  spherical  harmonics  of  e  which  will  now  be  proved 
and  which  permeates  many  discussions  to  follow. 

i  Kz 

The  plane  wave  e  is  a  solution  of  the  equation 

V2H*  +  =  0  (128) 

with  axial  symmetry.  This  equation  can  also  be  solved  in  spherical 
coordinates;  its  solution  is 


*  =  Pe(cos6)f£(r) 


029) 


■  JU4JULU1  tVMH  l<  \9V 


~V  •  >  -  y  vT.L,J.AiVWF»' 
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if  f  is  a  solution  of  the  equation,  then 

“7  ~  (r*2  HF>  +  ^  =  0  (130) 

r  dr 

where  P^Jcose)  is  the  nth  Legendre  coefficient.  That  is, 

fQ  =  eiKz  for  +  (K^  -  — ^ )  f&  for  z  =  0,  a  less  general  case. 

r 

Equation  130  can  be  solved  in  series;  there  are  two  solutions,  one 
beginning  with  rZ  and  the  other  with  r"®'"^ .  Let  us  denote  by  f  (r)  the 

JC 

solution  Equation  130  that  is  bounded  at  r  =  0.  Then,  except  for  an 
arbitrary  constant,  f0(r)  is  determined. 

If  A  are  arbitrary  constants,  e  A. P  (cose)f, (r)  is  a  solution 
z  j,=o  1  1  £ 

of  Equation  128  and  we  know  that  this  is  the  most  general  solution  of 

Equation  128  which  has  axial  symmetry  (i.e.,  does  not  involve  <j>)  and 

which  is  finite  at  the  origin.  [Note:  The  key  is  that  V2  ■*  (spherical 

Laplacian),  a  more  general  case,  with  the  restriction  that  f  is  bounded 

at  the  origin.  (Being  bounded  at  the  origin  implies  that  f  (r)  has  a 

power  series  expansion  such  that  the  pole  £^+-^  is  eliminated.] 

r2  ' 

i  Kz 

It  follows  that  e  can  be  expanded  in  the  form: 

eiKz  .  eiKrcos0  =  l  A  P  (cos9)f  (r) 

%•  0  *  1  (1 31 ) 

To  obtain  A^,  we  multiply  both  sides  by  Pa(cose)sine  and  integrate 
0  to  7r.  Putting  cose  =  t,  we  obtain 

1  •  ,  I  00 

/  elKrt  Pjl(t' )dt'  =  /  l  VW^P^t^dt*  (132) 

•*1  “1 

1  2 

Since  /  P£(t)P^(t: '  )dt'  =  ^f+T*  by  ortllo9onal'ity  of  Legendre 

functions,  we  have 


-1 


e1ICrtPA(t)dt 


2 

2Z+1 


VA(r) 


0  33) 


46 


.v\. y  - .>  ... . .. 


>>  y  .%  .v , 
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therefore 


1  JKrtP 


2I+TT  Wr)  *  I]  6  PA^dt  * 


2(i)£sin(Kr-  ^f) 


sin(Kr-^)  2A 

f£(r)  «  - ^ - —  and  ^  =  2(1  )£  or  «  (2*+l)(i)£ 


therefore 


e1Kz  *  l  (2t+l)(i)*  P  (cos9)f„(r) 


,  =  0  f0  '  ^  ■  J0‘Kr> 


f/  (r)s(2«F)1/2  Jt+|/2(Kr>  * 


We  can  look  at  the  general  solution  of  Equation  127 


(137) 


(138) 


039) 


V2iJ/  +  [K2  -  U(r)]ip  =  0 
2  2 

[NOTE:  V  =  and  is  not  necessarily  a  spherical  operator  (127) 
here]. 

As  before,  the  general  solution  of  Equation  127  having  axial  symmetry 


4>  =  l  AfP.(cos0)Lp(r) 


(140) 
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where  A  are  arbitrary  constants,  and  L?  is  any  solution  of 


(“V  ~  W-)  *  U2  -  u(r)  -  ^y^-)  1-0  =  0 

r  dr  rc  x 


(141) 


As  before.  Equation  141  has  two  independent  solutions,  one  finite 
at  the  origin  and  the  other  infinite. 


We  wish  to  choose  the  constants  A  so  that  Equation  140  shall 


represent  an  incident  wave  and  a  scattered  wave,  i.e.,  so  that  Equation 

iKr  PiKr 

140  shall  have  the  asymptotic  form  of  i{i  ~  e  +  - - f(e),  It  is 


necessary  that  our  wave  function  should  be  finite  everywhere,  L  must 


therefore  be  chosen  to  be  a  solution  of  Equation  141  that  is  finite  at 

G,(r) 

the  origin.  If  se  set  L£(r)  =  — - -  , 


dLfc  1  dG* 


r  dr‘7 


_  i  .  i  .  2G«. 

S7""7^  'ST-' 7S7 +  7" 


2dG*  ,  i dG*.  „  2h 

'S7"  T- 


Substituting  into  Equation  141  we  have 


1  (It  dL& 
7<2r-dT 


+  r2  +  (K2  -  U(r)  -  =  0 

dr  r 


1  dGt  2rG* 


I  I  _ *  _ *  ,  l-c  **  ,  I 

T  2r  '  7  SF-~r*  r  'TF*; 


2  (-2  ^  +  1  d  + 


+  ~T] 
r3 


+  (K2  -  U(r)  - 


*  ^  rp-’j 


’l 
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or 


d2G 


dr 


—  +  (K2  -  U(r)  -  G?  =  0 


(142) 


For  large  r,  the  last  two  terms  in  the  parenthesis  tend  to  be  zero, 

thus  we  should  expect  that  the  asymptotic  form  of  any  solution  G  would  be 
d2 

a  solution  of  -/  *  K2Gt  -  0  or  Gt  ■  A,cosKr  +  B,  sin  Kr  *  A  sin  (Kr  +  O 
dr 

where  A  and  e  are  constants. 


To  test  whether  this  is  so,  we  set 

G  =  U(r)e1Kr 

substituting  into  Equation  142  we  have 

«  u*elKr  +  UikelKr 
dr 


(143) 


or 


or 


*  U"elKr  +  2iKU'eiKr  -  K2UeiKr 
dr 


U"  +  2i KU '  -  Kjj  +  K2U  -  Uu  -  U  *  0 

u  r 


U"  +  2i  KU  *  -  Uu  -  A^+l)l  =  0 

rc 


(144) 


u"  +  ^  -  [u  +  U  =  o 


(145) 


For  large  r,  we  may  assume,  since  u  is  nearly  constant,  that 


d‘u  ^  KdU 

T7  Sr' 
dr 
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Neglecting  the  former  term,  we  can  integrate  Equation  145 

2rK  -  [u  +  *^^1]  U  =  Q 
ar 

2rk  .  [u  ♦  ilS$U]  <,r 

u  r 

Integration  gives 

> 

2iK  log  U  =  /r  [u(r)  +  dr  (146) 

r 

/.  u  =  exp  ik  /r  tu(r)  +  £^*+11]  dr  (147) 

r 

The  right-hand  side  of  Equation  145  tends  to  a  constant  for  large  r 
if  and  only  if  u(r)  tends  to  zero  faster  than  ^  as  r  tends  to  infinity. 
Thus,  for  fields  which  fall  to  zero  faster  than  the  Couloumb  field,  G, 
the  asymptotic  form  is:  G  *  A  sin(Kr+e). 

The  particular  solution  of  Equation  142  that  is  finite  at  the  origin 

Csin(Kr  -  1/2  JU+N 

will  then  have  the  form - ^ - -  where  C  is  an  arbitrary 

constant  and  N  is  a  constant  that  depends  on  K  and  u(r).  [NOTE:  the  term 
-  1  fU  is  added  so  that.  If  u(r)  is  zero,  N^  shall  be  zero.  The  reader 
may  feel  that  this  is  handwaving,  however,  this  form  shall  be  proven  in 
Appendix  A  through  perturbation  theory]. 

We  have  to  choose  the  constants  A  in  Equation  139.  If  we  subtract 
the  expansion  for  eiKz  for  the  Incident  wave,  we  obtain  the  expression 
for  the  scattered  wave.  We  shall  use  the  following  equations, 

OD 

^  *  l  A  l4(r)P  (cosG)  (148) 

i,=0  *  *  * 

s1n(Kr  -  1/2  Air  +  N£) 

K r 


L*(r)  = 


(149) 
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elKz  =  l  (2£+l)  1AP&(cosO)ft 


£=0 


*  *  eiKz  ♦  tm 


iKr 


(150) 

(151) 


together  with  the  fact  that  there  must  be  no  terms  proportional  to 
.-iKr 


— - —  in  the  asymptotic  expansion  of  the  scattered  waves  since  these  can 
be  shown  to  correspond  to  incident  waves.  Therefore,  looking  at  the  nth 
component : 

,iKrr 


.  e  p 

A  L  (r)  *  elKz  +  r  - ^os8)  =  eiKz  +  f(Q)e 


i  K.r 


(152) 


then 


JKTp 

Vt(r)  -  elKz  *  C*  - 


(153) 


or 


eiKr-l/2Zni  +  N  1  e-iKr+ i/2Zn  -  iN0  . 

jA£  - 2iKr - -  A£  - m? - -  (2^1  )l\} 


,iKr 


x  Pj(cos8)  =  Ct -p  ,t 


Mcos0) 


since  fff  =  sin(Kr  -  = 

*  RY  d 


ju.  Tr£i  -IKr  +  tt£ 

e 1  Nr  -  g  -g- 

2i  KT 


(154) 


substituting  this  into  Equation  154  we  have 


eiKr  -  If1+NAi 

li  KT 


-  (2£+l )  i 


£  e 


i  Kr  -  tt 


{ 


e“  1Kr  +  l§2--iN. 

A  - — = - -  -  (9o+i 

Xf 


rrtf — *-  <a+m: 


p,,(cos9) 
JTttr  }Pj(cos®) 


2Tkt, 


-iKr + 
e  c. 


(155) 


3  C£  ^P£(COS0) 
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JKr  -  nv’i 


•  iKr  +  i 


Matching  terms  of - 2 - an<*  - 2 


we  have 


J_  / 

2iK  \ 


IN,  ■  C  e 

A£e  1  -  (24+l)1AVe  x  *- 


r  -N0i 


2iK  1  £ 


-  (2A+1 )iA|*  0  from  Equation  157 


.  IN. 

A*  =  (24+1)1*  e  * 


Substituting  Equation  157'  into  Equation  15ttwe  have 


i*  (^)(e2iNl-l)e' 


-  irg.1 
2 


In  Equation  158  the  terms  i  and  e  can  be  combined  as 


TTj tJL  _ 

(e  2  )(e  2  )  -  1 


(||^)(e2i\l)  -  ct 


Therefore  the  final  equation  from  Equations  148  and  151  is 


<1 >=  l  (2£+l)i\|i(r)Pjl(cose)  since 


.iKr  C„ 


we  have  at  last 


f{0)  a  m  Z  (24+1 )(*  A-l)P4(cosG) 


(156) 


0  57) 


(158) 


059) 


(160) 
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This  last  expression  gives  us  the  expression  for  the  amplitude  of 
the  scattered  wave.  It  should  be  noted  that  f(o)  is  complex;  the 
scattered  intensity  1(e)  is  given  by  the  square  of  the  modulus  or: 


,2,21N* 


ft(9)f*(0)  = <2H+I)‘(e  1  )(e  1-l)pjj(cos0) 


since 


2iNp  -2iN?  2iN.  -21N 

(e  *-l)(e  -  1)  -  1  -  (e  +e  ;  +  1  -  2  -  2cos2N. 


iV0)!2"  “T  (2A+1  )2(1  -cosZN^) 
2K 


integrating  over  all  solid  angles 


Q=  l  2it  r|f.(0)|2sin0d6  =  J  (2-Z+]£‘-2  sin2N.P2 (cose)sinedQ 

H=0  JO  *  1=0  JQ  2K  Z  Z 


_  w 

2Tr  r- 

~1  L 

K  £=0 


l  (2il+l  )2sin^.  f  P2(cos9)sin9d9 
=0  n  £ 


since  the  integral 


f  P2(cos9)sin9d9  =  vie  have  Q  =  %  7  (2i,+l)sin2N 
JO  %  d!L+l  K  £=0 


[NOTE:  The  author  feels  that  for  a  complete  exposition  of  scattering 

theory  it  was  necessary  to  derive  the  steps  leading  to  the  scattering 

amplitude  Q.  This  section  will  eventually  lead  into  the  Generalized 

Theory  of  Scattering  which  is  the  crux  of  the  theoretical  foundation 

laid  for  the  understanding  of  the  computational  approach  to  the  solution 

of  the  carbon  complex.  However,  | f (e ) J  or  the  scattered  amplitude  is 

not  addressed  in  the  computer  programs  nor  is  it  used  explicitly  in  any 

2 

of  the  calculations.  If  | f (e ) |  were  used,  it  would  define  the  in  situ 
location  of  charge  densities  in  the  carbon  complex;  however,  for  multi- 
center  computations,  |f(e)|  would  be  untenably  complicated.] 
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1.  THE  BORN  SERIES  AND  T-MATRIX 

The  standard  problem  of  the  one-particle  in  quantum  theory  is  the 
effect  produced  upon  a  free  particle  by  a  local  potential  V(r^). 

It  has  been  taught  that  we  should  tackle  the  problem  by  looking  at 
solutions  of  the  corresponding  Schrodinger ' s  equation  with  various 
boundary  conditions.  However,  the  Green's  function  method  if  often  more 
powerful.  Since  it  has  already  been  discussed,  we  shall  continue  in 
this  respect. 

Looking  at  the  Schrodinger  equation 

H*  *  e*  and  O').  HQ<i>m(r)  = 

then  subtracting,  we  have 

(e-H)ty=  0 

Considered  as  an  operator,  the  Green's  function  satisfies 

(e-H)G  *  1 


(161) 

(162) 

(163) 


=>  G(e) 


1  _  1 
”  H-e  "  6  -  Hn  -  V 


(164) 


if  V  were  in  some  sense  "small",  this  could  be  expanded  in  an  infinite 
series 


^ i A HTCK .TOW .TPJtf JJJUgJRJBJJIJi WfflW  -V  U3W  j.-  .» 


1 


S.."VI 

„-'4 

F.’.V 


AFWAL-TR-83-11 30 


where 


1 


a  G-  +  G  VG  +  G  V  G  V  G  +  . . .  Where  GA  = - rr 

0  0  0  0  0  0  0  e  -  Hr 


letting 


G  >  G0  *  G0  V  G 


(166) 


since 


vs*«*  vv{vs>»}{vv«}«* 


(167) 


Equation  166  is  essentially  an  integral  equation  from  which  the 
Green  function  G  may  in  principle  be  derived. 


We  are  concerned  with  the  scattering  problem  in  which  the  potential 
V(r)  is  localized.  It  then  becomes  important  to  construct  Green's 
functions  that  refer  to  natural  casual  conditions,  with,  say,  a  state 
|i|>+>  that  looks  like  a  simple  forward  propagating  free-space  function  <|> 
at  a  large  distance  from  the  scatterer. 


These  functions  satisfy  the  integral  equation 

|ip+,  r>  *  |0,  r>  +  /sj(r,  r\  e)  V  (jr ' )  |/,  r’>  d3r' 


(168) 


This  equation  is  the  representation  in  real  space  of  the  Lippman- 
Schwinger  equation 


!♦*>■  l»>  +  rn^rrfv  l*+> 


(169) 


The  iterative  solution  of  Equation  169  is 


i*+>  ■  i*>  *  rro; v  !►  +  r-~}+  a » ,-tot  ’ » 


■  {i  +  gJ  v  +  gJ  v  G0  v  ♦ 


(170) 
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In  the  context  of  the  scattering  problem,  this  is  called  the  Born 
series.  The  term  of  nth  order  in  V  corresponds  to  successive  virtual 
scatterings  of  the  particle  by  the  scatterer.  To  bring  this  sort  of 
calculation  down  to  earth,  let  us  evaluate  the  first  term.  For  a  plane 
wave  of  momentum  K,  incident  on  a  localized  scatterer.  Equation  168 
becomes : 


4.  i  iKjr-r' |  , 

*  (r)  =  e  K  -  -T-  t  — vi — V(r')</>  (r)  d3r' 


(171) 


The  iterative  solution  of  Equation  171  requires  successive  integra¬ 
tions  over  all  space.  But  the  first  Born  approximation,  where  the  unper¬ 
turbed  plane  wave  is  put  in  place  of  the  true  scattered  wave  under  the 
integral  sign  can  be  evaluated.  In  scattering  problems  we  are  only 
interested  in  the  outgoing  wave  at  a  great  distance  r^  from  any  part  of  the 
region  where  the  scattering  potential  is  appreciable.  Asymptotically,  we  have 
i K | jr-r  *  |  iKr  -i(k/r)r*r' 


.  ,2 

r.  r '  r 

*  r  (1 - 2 - 2~)  ^  r  -  (p  •  r'  or 

r  r 


i  K |  r- r.'  I  iKr-K(r/r)*jr' 
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The  result  is  then: 


+(r)  *  .  _L  f  e!_!£] 

-  e  4it  J  jr-r' 


i  r  iK|r-r'|  ^ 

—  '  e - V(r').*1-1  d3r' 


ix-r  JKr  i  f  i (K— K* ) * r* '  . 

:0'HL_e - J-]v(r')*  “A' 


r  4ir 


dr1 


(172) 

(173) 


where  K'  is  a  vector  in  the  direction  of  r,  but  having  the  same  magnitude 
as  the  incident  wave- vector  K. 


The  term-by-term  summation  of  the  Born  series  is  a  very  inefficient 
way  of  calculating  transition  probabilities,  unless  the  convergence  is 
very  rapid.  This  is  not  always  the  case.  Almost  the  only  procedure  that 
works  is  to  separate  the  Schrodinger  equation  in  spherical  harmonics  and 
then  integrate  the  radial  function  outwards  to  evaluate  the  phase  shift 
N^(e)  for  the  partial  wave  of  angular  momentum  i  at  the  energy  e. 

The  exact  solution  of  Equation  171  at  large  distances  is 
i,  ji/«M  pi  Kr 

(r)  *  e  —  —  +  -  f(e)  (174) 

as  derived  earlier  and  where 

■,  00  21 N- 

f(9)  =  9Tk  l  (2t+l)(e  -1  )Pj  (cos9).  (175) 

£=0  * 

Comparison  with  Equation  173  shows  that  the  first  term  in  the  Born 
series  is  equivalent  to  making  the  approximation 

47rf(0)  <0K,  ! V|  (176) 

from  comparing  Equations  175  and  174. 

We  now  deliberately  invert  this  relation,  to  ask  the  following 
question:  what  operator,  t  says, gives  the  exact  scattering  amplitude 
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(Equation  176)  for  its  matrix  elements  between  the  free-particle  states 
|<t>K>  and  1 4>K , > ?  We  study  the  same  properties  of  this  transition  matrix 
or  t-matrix,  which  is  defined  as 

<0^i  1 1 1  =  -  4irf (9)  (177) 

In  other  words,  t  is  what  V  would  have  to  become  if  the  Born  formula 
(Equation  173)  were  to  yield  the  correct  transition  probability  or: 

<0K'  *1?  =  <0K'  ^  (178) 

or 

t\fiK>  =  V|/K>  079) 

Here  (<(>  >  represent  a  free-particle  wave  function,  which  acquiring 

lx 

an  outgoing  ?cattered  wave  when  modified  into  the  function  |^K>. 

The  significance  of  the  t-matrix  is  that  for  the  potential  scattering 
by  a  spherically  symmetrical  field,  the  angular  momentum  representation 
is  especially  favorable,  for  the  t-matrix  is  then  diagonal.  In  other 
words.  Equations  174  and  177  are  equivalent  to  writing  for  r,  r'  large 
as  t(r,  r',  e)  =  I  t^  (Kr)jJL(Kr')  Y^rjY^r)  where  j^  is  a  spherical 
Bessel  function  and  Y^r)  a  spherical  harmonic  for  the  direction  of  the 
vector  r.  The  diagonal  elements  of  t  in  this  representation  are  given 

V  s  ?tk  (exp  21  V'>  -  '>• 

The  t-matrix  will  again  appear  in  the  theory  of  generalized 
scattering  to  follow  -  so  this  brief  discussion  is  justified. 
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SECTION  VIII 

MULTIPLE  SCATTERING  METHOD  AS  RELEVANT  TO  THE  COMPUTATION 

The  multiple  scattering  method,  in  our  case,  deals  with  nonoverlapping 
spheres,  and  it  includes  not  only  these  spheres  and  the  interatomic  region 
between  them  but  also  an  outer  sphere,  surrounding  the  whole  molecule,  in 
which  the  wave  function  is  again  expanded  in  spherical  harmonics  and  radial 
functions. 

The  problem  we  take  up  is  conceptually  simple.  It  is  that  in  which 
the  potential  is  spherically  symmetrical  within  nonoverlapping  atomic 
spheres,  but  equal  to  zero  in  the  space  outside  these  spheres  and  out  to 
infinity. 

The  region  inside  the  spheres  is  designated  as  Region  I,  and  that 
outside  as  Region  II.  In  region  II  Schrodinger's  equation  is  the  ordinary 
wave  equation. 

V2U  *  -eU 
or 

£*  _[e . 

dr  tc 

where 

u  - 

We  are  dealing  only  with  negative  energies,  or  bound  states  of  the 

p 

molecule.  The  types  of  solution  we  wish  to  use  for  — ,  in  spherical 

m  \  "  _ 

coordinates  are  called  K^v  '  (Kr)  and  i£(Kr)  where  K  =  v-z.  These  functions 
have  forms  listed  in  Appendix  B. 

These  formulas  for  K  ^  (Kr)  Include  the  decreasing  exponential, 

X* 

exp(-Kr),  so  that  they  all  go  to  zero  at  infinity.  The  formulas  for  i^(Kr), 
being  regular  at  the  origin,  can  be  expanded  in  power  series  in  r. 


(180) 

(181) 
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In  Region  II,  which  extends  out  to  infinity,  we  can  build  up  the 
solution  to  Schr*ddinger's  equation  by  superimposing  functions  of  the  form 
^(Kr)Ym(e,<j>)  around  each  nucleus.  This  function  will  be  bound  to  go 
to  zero  at  infinity,  and  being  a  linear  combination  of  solutions  of 
Equation  180,  all  for  the  same  energy,  it  will  be  a  general  solution  of 
that  equation,  going  to  zero  at  Infinity.  That  is,  as  the  general  solution 
in  Region  II,  we  may  write 

*n  =  082) 

where  B  denotes  the  Bth  atom,  |rD|  is  the  distance  from  the  Bth  nucleus 
to  the  point  where  the  wave  function  is  being  computed,  and  Y,mCrn)  means 
the  spherical  harmonic  for  the  angle  of  the  vector  fg,  using  the  Bth 
nucleus  as  origin.  We  must  now  apply  the  boundary  conditions  to  make 
this  function  and  the  solution  in  Region  I  continuous  with  continuous 
derivatives  over  all  the  spherical  surfaces  bounding  the  atomic  spheres. 

U(e,r,9,0)  =  (183) 

for  the  equation 

V2U  =  -(e  +  V) U.  (184) 

Thus,  U  is  a  solution  of  the  Schrodinger 's  equation  regular  at  the 

origin,  containing  an  infinite  number  of  arbitrary  constants,  the  C  's. 

am 

The  type  of  expansion  of  the  wave  function  (Equation  182)  is  one 
possible  way  of  describing  it.  An  alternative  form  of  expansion  is, 
however,  more  convenient  to  use  in  treating  the  continuity  of  the  wave 
function  over  the  surface  of  the  spheres.  Suppose  we  wish  to  expand 
within  a  cell  containing  the  particular  atom  a.  There  will  then  be  two 
types  of  terms  in  Equation  182;  those  for  which  b  =  a,  the  functions 
representing  the  tail  of  the  atom  a^,  and  those  for  which  b  f  a,  repre¬ 
senting  all  the  other  atoms.  We  shall  expand  these  tails  of  atoms  b, 
within  the  cell  a,  in  terms  of  r,  and  spherical  harmonics  of  the  angle 

a 

r  .  We  employ  a  result  of  a  general  theorem  by  which  we  can  expand  a 
a  #  1  \ 

function  of  the  form  K„  (KrD)Y„  (rD)  as  a  linear  combination  of  functions 

%  □  dm  b 

i  ,(Kr  )Yn,  (r  )  around  the  center  a. 

3  x,  m  a 
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Since  the  function  i  (Kr  )  can  be  expanded  in  power  series  as 

=  1.3-5-^211+1)  }  +  2(2£+3)  +  2-4- ^2A+3)(2£+l )  +  **]  ( 

and  since  spherical  harmonics  ,m i (r  )  can  also  be  expanded  in  power 
series  in  the  x,y,z  coordinates  around  nucleus  a^  the  functions 
i  .{Kr  )Y  (r  )  can  also  be  expanded  in  x,  y  and  z. 

£  d  £ITI  a 

This  theorm  is  (Reference  3) 


(185) 


+  N  ft  ‘m1  •  ) 


^KraB^YL,m'-m^raB^  V  *Kra^Vm' *ra* 


(186) 


raB  =  ra  -  rB 


(187) 


IL(tm;£n,')  =  /YL>m.m,  (r^rjY^,  (r)dfi 


(188) 


where  the  integration  is  over  solid  angle. 

The  integrals  over  a  product  of  three  spherical  harmonics  are  the 
same  ones  met  in  the  definition  of  the  Clebsch-Gordon  coefficients  in 
the  theory  of  atomic  spectra.  (They  are  treated  in  more  detail  in  the 
section  on  generalized  scattered-wave  theory). 

Let  us  write  the  whole  function  in  this  form.  We  then  sum  over 
all  nuclei  b  except  the  ath  and  include  separately  the  term  for  nucleus  a. 


We  have 


*11  =  E(b,£’m^ALnK£1^KrB^YJUn^'B^ 


*  E[ALKi1)(Kra)V?a)  +  2  A£m  4it(-1}* 

£m  m  *  a  urn  a  b/aA'm'L  m 

4 

x  ILU'">';to)K[1)(KraB)Y*>m,.m(faB)1ll,(KrJ)Yl,|I1,(?a)] 


(189) 
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w  \TVO»^  V*  w.->» 


mi  'A  ',*  ■'  -  '.iff-  'A'AITA  ’■■  \.\  A  . 


In  the  second  term,  we  interchange  the  labels  t,m  with  z',m'.  Then  we 
have 


*»  ■  J v- 


(?a)  KrfH* 


+  £  A?»mi  I,  (ZmjA'm' ) kA 

B^a.Z'm'L  1  m  L  L 


i*  0  'm*  \  ) 


x  <*<»'.>]■  *  W 


(raB)AVitT  VKra* 


where 


Zm.S/m1  aB1 


aB 


x  YL,m-m' 


We  can  rewrite  Equation  189  as 


*11  =  ?„  V^Il'11"0 

zm 


where 


♦  jjUm) 


is  the  bracket  multiplying  Y£m(ra)  in  Equation  190. 
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(190) 


(191) 


(192) 


■JJURVJMJl  V  U 
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We  must  now  produce  continuity  of  the  wave  function  and  its  first 
derivative  immediately  inside  the  sphere  surround  the  ath  nucleus,  where 


the  wave  function  is  a  linear  combination  of  functions  R„(Kr  )Y.  (r  ), 

t  a  tm  a 


and  immediately  outside  the  sphere,  where  the  function  is  given  by 
Equation  189.  For  continuity  of  function  and  derivative  at  the  surface 
of  the  sphere,  we  make  the  logarithmic  derivative  of  the  term  for  each  t,m 
continuous  as  well  as  the  function  ip  itself  and  its  normal  derivative. 

That  is,  at  the  radius  of  the  ath  sphere: 


1 


#n(£m) 


VjjlAm)  dr. 


dr 


(193) 


gives 


^rI(4m)  dR? 

R£(Kra)  -  d/jj (Am)  ^ -  (Kra)  =  0 

3 


if: 


(194) 


If  we  use  bracket  notation  for  the  Wronskian,  defined  as 

Rfvi  m*± 


[A(x),B(x)3  -  A(X)  B(x) 


(195) 


we  can  rewrite  Equation  194  as 

-.a 


CR®(Kra),  ^XI(Am)]  *  0 


(196) 


When  we  evaluate  the  Wronskian  (Equation  196),  we  will  have  at  first: 

r(D 


— A*m  l  (-I^VgJ*  ,  , 

£fl|  Qra  3  R*aO»  m'  m 


Hr 


Bj*a£'  ,m' 


(faB’e)  A£\m'  af^ 


di 


(197) 
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dK/1}  .  di,  ,  m  dR* 

This  Involves  terms  like  —  •  R*  -  —•  R*  -  '  ~ 

a  a  a 

dRn 

+  \(^ra)  (neglecting  other  coefficients  not  a  function  of  rg  for 
a 

.  f^L11  a  m  dRil  rfdiL  dRtl 

simplicity.)  This  is  R4  -  K^*  jr-]  -  R,  VKra>dH 

a  a  a  a 

or  ^R®,  K^j  -  \(Kra)j.  Putting  this  into  the  total  Wron- 


skian,  we  then  have: 


[Rj,  K<’)(Kra)] 


B,.Jv  <Vd> 


(198) 


x  (-1)il'A|.n).  CRj(Kra).  1t(Kra)]  -  0 


Dividing  both  terms  of  Equation  198  by  [R®  (Kr*a),  1 & ( kf  ) ] ,  we  have 


(a  HRg(K.r-a),41)(Kra)3 

*■  [RJ(Kra>-  V(Kra> 


(199) 


-  I 

B^a.A'.m' 


K 


-lraB 


;e)(-l) 


where  we  have  such  an  equation  for  each  a  (i.e.,  for  each  atom),  and 

for  each  t,m,  for  which  the  coefficients  A®^  are  of  appreciable  size. 

Equation  199  gives  us  the  coefficients  A^m  of  the  wavelets  in  Region  II. 

We  need  in  addition  the  coefficients  C  of  Equation  183  for  the  expansion 

s,m 

inside  the  atomic  spheres.  Let  us  derive  the  relationship  between  the 
A's  and  C's. 
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We  use  Equation  198  to  express  the  summation  over  b,  e,  m'  which 
appears  in  Equation  190  in  terms  of  A®m.  We  demand  that  the  summations 
of  Equations  190  and  183  representing  the  solutions  inside  and  outside 
the  ath  sphere,  should  be  equal  at  the  radius  r  of  the  sphere. 

d 


This  demands  that  the  coefficients  of  the  spherical  harmonics  Y£m(0»4>) 
in  both  summations  should  be  equal  for  each  z  and  m  as  well  as  for  each 
sphere.  Then  using  Equation  199  we  have 


CR>'ra>  *  ALKinCKra)  - 


[^(a-r^.K^CKr)] 

[Rj(«.ra).1t(Kra)] 


VKra> 


since 


(200) 


dj(P 

i£(Kr)  d(Kr) 


(Kr)  -  K^Kr) 


diA(kr) 
d  (Kr)' 


We  can  reduce  Equation  200  as  follows: 

cLRa(‘"ra>  ■' ALKi1>tlCra){ 


Rad1, 

W~ 


[RMA(Kra)3 


{ 


dR 


.a*!0  .(!)_, 

K«.  dra  *  h  dR 

a  o 


■}  VKra> 


CR'.i.tKra)] 


5l 

'tdrj 


(201) 
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0$*W3 


dSt 

♦  drj  i 


K(1)i  -  Kll)V 

l  h  h  >1 


i} 


CR;,ittKra)3 


a'  a' 


i+l 


-  A 


,»  (-D, 

- *  c!„R!(e,rJ 


or 


al  =  cL  ^4-  crJ.  yKra)] 


£m  Jim 


(-i  r 


which  is  the  desired  relation  between  the  A's  and  the  C's, 


(202) 


The  above  equations,  and  in  particular  Equation  199,  constitute  the 
fundamental  equations  of  the  multiple  scattering  cluster  method.  The 

g 

coefficients  A^m  determined  by  these  equations  are  the  coefficients  of 
the  expansion  of  the  wave  function  in  the  interatomic  Region  II,  in  terms 
of  the  scattered  wavelets  emitted  from  each  atom.  As  we  have  seen,  for 
negative  energy  these  wavelets,  instead  of  being  sinusoidal  as  they  would 
be  for  positive  energy,  are  falling  off  exponentially  as  we  go  away  from 
the  Bth  nucleus.  We  have  as  many  simultaneous  linear  equations  as  we 
have  wavelets.  In  many  cases  this  means  that  only  small  value  of  i  need 
be  used  because  the  expansion  in  scattered  wavelets  is  rapidly  convergent. 
In  general  the  equations  will  have  no  solutions,  and  it  is  only  if  the 
determinant  of  the  coefficients  of  Equation  191  is  zero  that  we  get 
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solutions.  Since  the  energy  e  is  entering  into  the  coefficients,  we  have 
solutions  only  for  certain  energy  eigenvalues  of  the  problem,  which  is  the 
method  by  which  the  energy  levels  are  found. 

For  each  eigenvalue.  Equation  191  can  be  solved  for  the  coefficients 
A^m  and  hence  for  the  wave  function  in  Region  II.  In  turn,  from  Equation 
202,  we  can  find  the  wave  function  inside  each  sphere. 


tttttjttttttttt:  j.  rj  ?.i 
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SECTION  IX 

GENERALIZED  SCATTERED-WAVE  APPROACH  TO  MOLECULAR-ORBITAL  THEORY* 

1 .  INTRODUCTION 

This  approach  is  an  adoption  of  a  rapidly  convergent  partial-wave 
representation  of  the  multi  center  wave  functions,  the  use  of  multiple 
scattering  formalism,  and  the  separation  of  the  principal  equations  into 
"structural"  and  "atomic"  quantities.  These  features  distinguish  this 
technique  from  more  conventional  LCAO  type  molecular-orbital  methods. 

This  Hamiltonian  is  based  on  the  partitioning  of  molecule  space  into 
regions  of  nonoverlapping  spherically  averaged  potential,  including  a 
local  statistical  approximation  to  the  exchange  potential.  There  are 
no  multicenter  integrals  to  compute  or  to  estimate  as  there  are  in  LCAO 
methods.  The  scheme  is  comparatively  easy  to  implement  computationally 
and  is  practicable  for  systems  of  considerable  stereochemical  complexity. 


2.  SOLVING  THE  INTEGRAL  WAVE  EQUATION 

Consider  an  N-center  molecule  having  atoms  located  at  the  sites 
R  (a  =  1 ,2,. . . ,N),  with  respect  to  some  arbitrary  origin.  In  the 

r  U 

scattered-wave  approach  to  molecular-orbital  theory,  we  are  required  to 
solve  the  integral  Schrodinger  equation:  (See  Appendix  C) 

tj>(r)  a  /dr'GQ(r,r' ;e)  V(r' )  Mr* )  (203) 

for  the  stationary  states  i|»(r)  and  energies  e  of  single  electron  bound  to 

the  local  potential  field  V(r)  of  the  molecule.  V  (r)  includes  the  Hartree- 

Fock-Slater  or  X  statistical  approximation  to  exchange,  along  with  the 
a 

average  Coulomb  potential.  We  assume  that  the  molecular  potential  at  an 
arbitrary  point  £  can  be  represented  as  a  superposition  of  overlapping 
local  potentials  centered  on  the  atomic  sites,  namely. 


V(r) 


N  N 

l  V(r-Rg)  £  l  V(rJ 
=n  8=0  ^ 


e=o 


(204) 


♦This  chapter  is,  in  large  part,  a  verbatim  synopsis  of  Reference  11  by 
K.  H.  Johnson.  This  author  is  gratefully  acknowledged. 


69 


\  "„  ".  >  *'t 


AFWAL-TR-83-1 1 30 


Gq(t,  r';  r.)  i s  the  "free-space"  single-particle  Green's  function 
satisfying  the  inhomogeneous  wave  equation. 

(V2  +  e)  GQ(r,  r';  e)  -  <5(r  -  r')!  (205) 

For  a  polyatomic  molecule  with  overlapping  potentials  of  Equation 
204,  the  integral  (Equation  203)  cannot  be  easily  solved.  However,  for 
the  special  case  of  nonoverlapping  spherically  averaged  potentials,  we 
can  use  Green's  theorem  to  reduce  the  volume  integral  in  Equation  203  to 
a  sum  of  surface  integrals  over  spheres,  leading  to  an  exactly  solvable 
model . 


In  order  to  develop  a  more  flexible  class  of  model  potentials  for 
which  Equation  203  is  still  approximately  solvable,  let  us  consider  the 
general  problem  of  determining  solutions  of  the  integral  equation: 

\p(r)  =  /dr*  Jdr"  GQ(r,  r";  e)  U(e;  r",  r')^(r')  (206) 

for  an  energy-dependent,  non-local  potential  operator  U(e;  r,  rl ). 
Equation  206  is  solvable  provided  we  expand  the  potential  operator  in  the 
separable  multi  center  partial -wave  representation, 

N  „  * 

U(e,  r,  r')  *  ^  £  U[(e)  ^(r^ )<*>L (jr^)  (207) 

where 

\(lg)  =  4>L(r6)  YL(r$>  (208) 

here  <i>  (r  }  has  been  expanded  in  terms  of  the  product  of  spherical 

L  “p 

harmonics  Y^(v,$)  and  radial  amplitude  <^(rg). 

Also,  in  Equations  207  and  208  L  =  U,m)  is  the  partial-wave 
(angular  momentum)  index  and  U®(e)  represents  certain  energy-dependent 
coefficients;  U®(e)  and  $L(rg)  are  chosen  t0  simulate  (exactly  or 
approximately)  the  originally  defined  molecular  potential  (Equation  204) 
or  more  general  potential  configurations.  Both  may  be  complex,  in 
general,  to  allow  for  the  possible  use  of  "optical"  or  "self-energy" 
potentials. 
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\'  .  \  v  ■"'Tf  7*  ’  \-  ^c^,3!'T^T^jT^r5T^*y 


Solutions  of  the  integral  equation  can  be  generated  by  substituting 
Equation  207  into  Equation  206  and  expanding  the  Green's  function 

*  £)  anc*  wavefunction  ip (ir )  in  the  partial-wave  representation. 


A  convenient  integral  representation  of  the  free-space  Green's 
function  is  (see  Appendix  D): 

tyr,  i‘ ;  «)  *  ,-W  da  e»P^-  (y‘)  (jo,) 

(2tt)  e  -  q 

Because  GQ  depends  on  the  difference  between  "field"  and  "source" 
vectors,  we  can  choose  to  write  these  vectors  with  respect  to  the  atomic 
sites  R^,  namely 

-  expCri'(c-^-ILJ 
u  a  -B  /o-'~  e  -  q 


—a 


(2tt)' 


-Q 
T 


in  which 


R  0  =  R-  -  R 
-aB  -a 


Pictorial ly,  we  have 


Figure  2.  Diagram  Showing  Two  Potential  Points  r  and 
Two  Atomic  Sites  R 


(210) 


iQU 


l 


sa 


P 
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The  vector  R  +r  -  [R„  +  r„]  =  r  -  r '  -  R  .  represents  the  dis- 
—a  “a  -iJ  ~ P  -a  S  -otP 


placement  between  two  arbitrary  potential  points  with  respect  to  two 
arbitrary  atomic  sites.  We  shall  use  Equation  210  in  Equation  206  for 
those  terms  in  the  sum  over  atomic  index  6  for  which  &  /  a.  For  the  case 
where  6  =  a,  we  shall  use  the  alternative  representation  of  Gq,  given  by 
(Appendix  E) 


GQ(J1.  r';  e)  = 


J_  exp(-K|£.-r.' 

4tt  i  _  _i  i 


| r-r' | 


(211) 


where  for  bound  states  of  negative  energy  K  h  The  appropriate 

partial-wave  expansions  of  Equation  211  are 


G0(H*  I';  6)  =  -K  Z  (-l)\(Kr)  T^(Kr')  for  r  >  r* 


(212) 


Gq(£»  r';  e)  =  -K  I  (-l)£IL(Kr)  K^(Kr' ); 


r'  >  r  (213) 


In  these  equations 


IL(Kr)  =  iz(Kr)YL(r) 


(214) 


KL(Kr)  =  K[1}(Kr)YL(r) 


(215) 


where  i  and  kJ  J  are  modified  spherical  Bessel  and  Hankel  functions, 

X  X 


respectively,  defined  by 


iA(x)  =  1  j^O’x) 


(216) 


K«0)(x)  =  "i”^hP  }(iX) 


(217) 


«  ,  •  *•  .  •  ■-  ,>  • 


’  *  * '  .  •  *  ■  •  1 
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4>(r  +  R  )  =  EC?  R?  (r  )  y,  (r  ) 
-o  -a'  ^  L  L  '  cr  L'-a' 


(219) 


where  the  R“(r  )  are  unknown  radial  amplitudes  for  the  one-electror 
La 

orbitals  and  Y^r^)  are  spherical  harmonics. 

We  next  expand  the  plane  waves  in  terms  of  radial  functions  (Bessel 
functions)  and  spherical  harmonics  so  that  the  angular  functions^  are  all 
expanded  in  terms  of  the  same  basis  functions,  namely,  the  spherical 
harmonics.  Since  it  has  already  been  shown  that  plane  waves  can  be 
expanded  in  terms  of  analytic  Bessel  functions  explictly  as  well,  this 
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presents  no  problem;  however,  one  must  bear  in  mind  that  the  wave  function 

tj/(r  +  R  )  is  not  an  explicit  analytic  function  of  a  radial  coordinate, 

- a  ~a  - 

such  as  a  Bessel  function,  but  rather  is  a  more  complicated  radially 
dependent  function  to  be  determined. 


Expanding  the  plane  waves, 

exp(ia  'r,)  =  z  i\(qra)  Y*(a)  Y^) 

exp(-i£  •  r£)  =  4-rr  Z  i £ (qr * & )  YL(g.)  YL(r'e) 


(220) 

(221) 


and  using  Equations  208,  214,  215,  and  219  through  221  in  Equation 
218  we  have: 

\p  (r)  =  ip(r  +  R  )  s  E  C?  f**(r  )  Y.  (r  ) 

Y'— '  -a.  -a  ^  L  L  a  L  -a' 

(KrVY(^)  [.  *lK> 

\ _ l 

x  W  *  t<Kra>  W  JrfrjK  4V<K> 


[.  V(ra>  VL'<^>}  UL-<*>  IK  *L'(ra>  VL'<^> 


l <?  KW 


(222) 


m  •  4,  i  i\(qrX(9)  w  1  exp('yB> 


.  jewWiWiw  rrr^rr 
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x  f.  /dr*  V<re>  Yi^>  le  cl?  Rf’<re>  V<I$> 


From  the  orthonormality  of  spherical  harmonics,  we  see  immediately 
that  the  functions  under  which  brackets i— J  appear  can  be  integrated  over 
solid  angle  dfi  as  /dfl(£)  Y^(r)  Y^.(r)  =  6^1,  remembering  of  course  that 

K  C  fo  fa  d9“  sin9“  ir“r° 

Integrating  with  respect  to  the  spherical  harmonics,  we  have 


*(r)  =  t  Cj  R»(ro)  Yjr,)  .  -K  l  (-1  >‘{k<'  >(Kra>  Y^) 

X  Lr  dra  r ‘I  1*.(Kra>  *L(ri;)  +  W 


r’<r 
a  a 


X  lr.>r  dri  ra2  Ki1>(kri>  *L(rcl)}  «?(*)  C  K  rf 

a  a  u 

1  n  *  N  exp(-iq*R  ) 

*  toFM-4,Lt1  W  Va)  I  -  -r 


x  Y.  (r  ) 
L'-cr 


C  dr8  r82  J.  V<a>  V‘pe>  uL’(e>/dr-B  i 


X  4>L,(rJ)  C®,  R®t(rJ) 
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Looking  at  the  Lth  component  and  dividing  YL(r^)  out  from  both  sides  of 
this  equation,  we  have: 

c?  f«r.)  -  {t(l)(K'«>  f“  ri2  dri 

ra  a 

Wf  r'2  dr'  K^IKr')  +L<r;)}  C^(e) 


+  1, 


rct 
r>r 
a  a 


f0  r;2  <(•;)  $  *  f  «*/*  fa 


N 

x  l 

expf-ia-^g 

e-qZ 

Bfa 

x  C®, 

uL-(e>  C r; 

This  rather  complicated  looking  expression  can  be  simplified,  and 
the  terms  can  be  given  a  physical  interpretation.  First,  we  expand  the 
plane  wave, 

*’ip(-i*w  *  *"  l„  M)  \«(a)  (223) 

Inserting  this  into  the  corresponding  part  of  Equation  222,  we  have 

N 


1 ft a  V«r«>  vL(a>  i  4*  h  r" 

i  — i  PfCt  L  i — * 


X  ^..0^3) 


l  r*'  Cre2dr6^S'>\'‘re>  VL^! 


(a) 


■> 


/<«(a).YL„(9)  v^a)  yl,(i)  =  p-  (l-.l') 
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N  These  are  the  Gaunt  integrals,  which  are  non-zero  only  for 
1 4-2,'  |  <  i"  <  l  +  2'  and  i  +  2,'  +  2."  ■  even  integer. 

The  Gaunt  integrals  can  be  expressed  as  products  of  Clebsch- 
Gordon  coefficients,  and  tabulations  of  them  are  readily  available. 
Using  the  notation 


<?  <£<«>  f0  #r«*L<ra,l?(r«) 

(224) 

rl  “Vu'^a*  »L(rc,) 

(225) 

K(.,)W  £•  rf  K  ii(Kr.}  ^(r.) 

(226) 

(227) 

=  4u1  q  1  ^  k"^L;L‘)j’2:,(qRaB)YL"^e  (228) 
we  arrive  at  the  equation 

CLf?,‘ra)  =  +  9L(ra;e)  VKra)] 

+  %l  l*  #L(‘>  Fo  ,  dl>V‘>ra><*-«2>'1 

(229) 


m 
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Since  all  the  quantities  of  Equations  223  through  228  are  known,  the 
radial  functions  C“r“  (and  hence  the  effective  single-particle  wave- 
functions)  can  be  determined,  once  the  coefficients  A“  (Equation  224) 
are  calculated.  To  determine  the  A“  and  the  single-particle  energies, 
we  multiply  Equation  229  by  r^^r^)  and  integrate  over  all  r  .  Using 
Equation  225  and  collecting  terms,  we  obtain 


j,  [,  «cB«LL'  C^Ce)"1  -  |^rad'’c?L(rc,)fL(rc.ie)Ki1>(l<ra) 

-  f0  r»drA(ra)9L(rc,18)1l(Krcl)] 

■  <1-V  I  f  q  <Jn*aL<<))(e-<i2)‘,G“f.(^eiqML-(‘i)  **<«)  *  0 


(230) 


This  is  a  set  of  linear,  homogeneous  equations  in  the  coefficients  Aj\ 
To  give  some  physical  significance  to  the  matrix  elements  in  Equation 
230  we  divide  each  row  a  by  the  quantity 


5  C  #raVKrcXl’'a) 


(231) 


and  each  column  e  by  the  quantity 


E  j”  rs<Vl'(KrS)V(,'fS) 


(232) 


SC 

i: 


Then  we  can  write  Equation  230  in  the  form 


♦V‘«>  j,  l  -£.<«*«>}  *?.(•) 

x  <J)^(e)  =  0 


(233) 


araHwiiUBUK  mw  j.  a  abajjvafjiu  ^  ».t»  w«  gvir-  a  v*  w»  ^ 
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in  which  we  have  introduced  the  notation 

t|_(e)  5  (234) 

_ »n^(e)4f(«) _ 

[UL(e,r1  '  £  r2<lra{(ra)l:ft(ra;e)Kr)(Kra)+9L(ra;e,V(Kra)] 

1  (,‘ V  I  l>“L(e>*L,(6):r' 

(235) 

X  j“  q  <l<»“*(q)(e-q2)"'(^®,(^jaiq)^,(q) 

The  bound-state,  single-particle  energies  correspond  to  zeros  of  the 
secular  determinant  (Appendix  H). 

I4a6  •LL.ttJWr'  -^f.O^S.)'  (236) 

One  can  solve  for  the  coefficients  A^(e)  corresponding  to  these  energies 
via  Equation  233  and  then  for  the  effective  single-particle  orbitals  by 
Equation  229. 
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SECTION  X 

DISCUSSION  OF  THE  COMPUTER  PROGRAMS 

1 .  HSRHF  ATOMIC  PROGRAM 

The  program  of  F.  Herman  and  S.  Skillman  (Reference  5)  for  the  self- 
consistent-  field  Hartree-Fock-Slater  calculations  on  atoms,  has  been 
modified  to  allow  for  variable  alpha  in  the  X-alpha  approximation  to 
exchange  correlation. 

The  object  of  this  program  is  to  obtain  a  self-consistent  solution 
of  non-relativistic  Hartree-Fock-Slater  equations  for  a  free  atom  or  ion, 
given  the  atomic  number  Z,  the  electronic  configuration  (which  is 
specified  by  the  orbital  configuration  numbers)  and  the  self-consistency 
criterion  3  (SCF). 

The  input  information  includes  a  trial  or  starting  potential,  a  set 
of  trial  energy  eigenvalues,  as  well  as  other  information.  The  solution 
includes  the  self-consistent  potential  expressed  in  the  form  [rV(r)]; 
the  energy  eigenvalues  En  and  the  normalized  radial  wave  functions  PnX(r) 
at  all  (or  some)  of  the  441  points  of  the  integration  mesh  and  the  values 
of  PnX(r)/rA+l  ar  r  =  0.  This  information  is  expressed  in  floating  point 
form  to  eight  significant  figures. 

Some  of  the  data  required  includes  Schwarz's  alpha  parameter.  Green's 
parameter  for  generating  the  potential,  and  trial  energy  values.  The 
latter  are  obtained  from  Reference  5. 

The  calculated  free-atom  eigenvalues  and  charge  densities  are  input 
to  the  potential-superposition-program  MOLPOT  (Molecular  Potential). 

2.  MOLPOT 

MOLPOT  is  used  to  construct  the  initial  molecular  potential  for 
X-alpha  standing  wave  calculations,  by  superposition  of  free-atom 
potentials  and  spherical  averaging  of  the  superposition  in  each  atomic 
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and  the  extramolecular  region.  Free-atom  potentials  are  calculated  from 
the  input  HSRHF  charge  densities,  the  Coulomb  term  being  obtained  from 
Poisson's  equation  and  the  exchange  term  from  Slater's  X-alpha  approxi¬ 
mation. 


Various  calculations  are  required  as  input  data  for  M0LP0T.  For 
instance,  the  a  parameter  for  the  interatomic  region  is  needed.  This  is 
usually  a  weighted  average  of  the  atomic  o's;  possible  choices  for  the 
weights  are  the  atomic  sphere  volumes  or  numbers  of  electrons  per  atom. 

The  HSRHF  charge  densities,  one  for  each  unique  atom,  are  required; 
coordinates  of  each  atom  with  respect  to  a  symmetry  center  or  other 
geometrical  origin  must  also  be  provided.  This  calculation  can  be  involved, 
depending  on  the  complexity  of  the  molecule.  At  times,  a  physical  model 
of  the  complex  can  facilitate  geometrical  insight  into  the  relative 
positions  of  atoms,  nearest  neighbors,  etc. 


MOLPOT  provides  the  free-atom  potential,  charge  density,  and  integrated 


charge  p  -  f  4irr  p(r)dr,  printed  as  a  function  of  r  for  each  unique 


atom  prior  to  the  potential  superposition.  After  superposition,  the 


quantities  r,  V(r),  rV  (r),  V  p(r),  and  p  are  tabulated. 

c  r0  r0 


3.  HARMONY 


HARMONY  calculates  the  spherical  harmonic  basis  functions  for  a  given 
symmetry,  range  of  values,  and  arrangement  of  atoms,  by  a  projection 
operator  approach. 


Representations  of  a  given  molecular  point  group  including  inversion 
must  be  included  in  the  input.  The  basis  functions  which  are  generated 
are  then  orthogonal i zed.  This  output  is  punched  in  the  proper  format  for 
input  to  the  Energy  and  SCF  programs  (see  below).  The  calculations  are 
made  for  each  representation  for  each  group  of  equivalent  atoms  in  the 
molecule  which  have  different  angular  relations  to  each  other. 
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4.  ENERGY 

The  program,  ENERGY,  is  the  non-self  consistent  version  of  the  main 
X-aSW  program.  The  program  searches  for  eigenvalues  within  specified 
energy  ranges  for  an  input  molecular  potential  and  symmetry.  It  is  most 
commonly  employed  to  obtain  initial  energy  levels  to  start  the  SCF 
process,  and  to  convert  the  long- form  molecular  potential  from  M0LP0T 
to  the  short  form  used  as  input  to  the  SCF  program. 

It  is  used  to  search  for  excited  states  at  self-consistency,  or  to 
obtain  radial  wave  functions  and  their  coefficients  for  input  to  the 
contour-mapping  program. 

Input  data  include  radii  of  atomic  spheres  (including  the  outer 
sphere),  the  basis  functions  for  a  particular  representation  (from 
HARMONY),  and  the  energy  range  of  valence  levels. 

The  output  from  ENERGY  consists  of  (1)  the  input  potential  in  each 
region,  always  converted  to  the  short  form  in  print,  and  if  desired  in 
punch,  and  (2)  for  each  representation,  the  input  symmetry  and  a  table 
of  energy  values,  determinant  values,  and  Ramsauer  factors  in  the 
specified  energy  ranges,  with  indicated  determinant  zeros. 

The  Ramsauer  factor,  always  +1  or  -1,  is  included  since  the  program 
will  occasionally  Indicate  a  zero  due  to  a  change  in  its  sign,  when  in 
fact  the  determinant  has  not  changed  sign,  the  latter  being  the  criterion 
for  a  true  zero. 

5.  SCF 

The  self-consistent  version  of  the  main  X-alpha  SW  program  searches 
for  eigenvalues  at  E^  +  aE^  ,  where  E^  and  aE..  are  specified  for  each 
state  i,  and  n  takes  on  the  values  1,  2,  3,  etc.  until  a  determinant 
zero  is  found.  It  calculates  normalized  probability  densities  in 

the  various  regions  for  each  state  i_.  After  all  states  have  been 
calculated,  occupation  numbers  are  assigned  according  to  input  criteria, 
and  a  new  many-electron  charge  density  and  potential,  as  well  as  total 
and  kinetic  energies,  are  obtained.  The  new  potential  is  then  normally 
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used  as  input  for  another  iteration  until  the  eigenvalues  converge  to 
+  0.001  Rydbergs. 

Input  data  include  energy  eigenvalues  of  each  shell  for  each  atomic 
representation,  the  HARMONY  output  deck,  and  the  potential  deck  from 
MOLPOT. 

Printed  output  from  SCF  consists  of  the  input  potential  and  symmetry, 
followed  by  the  eigenvalue  determination,  normalization  factor,  and 
normalized  probability  densities  in  the  various  region  for  each  state. 

At  the  end  of  each  iteration  the  eigenvalues  and  probabilities  are 
summarized  along  with  the  total  and  kinetic  energies,  the  various 
components  of  the  total  potential  energy,  the  total  charge  in  each 
region,  the  largest  fractional  difference  between  the  old  and  new 
potential,  and  the  new  potential  itself. 


r?«  . 
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SECTION  XI 
THE  CALCULATION 


1 .  STATEMENT  OF  THE  PROBLEM 

The  cluster  complex  to  be  calculated  is  CgHg.  The  energy  eigenvalues 
shall  be  calculated  and  compared  with  those  of  CgH.^  obtained  by  Watkins 
and  Messmer  (Reference  28).  This  intercomparison  will  determine  whether 
the  lowered  symmetry  of  CgHg  (Cgv  symmetry)  shifts  the  eigenvalues  of  the 
more  symmetrical  cluster  CgH.^  (T^  symmetry).  Physically,  if  the  eigen¬ 
values  are  shifted,  it  may  happen  that  the  bulk  states  at  the  top  of  the 
valence  band  of  diamond  will  split  away,  forming  a  surface  state  within 
the  gap.  This  state,  if  it  exists,  must  be  associated  with  the  three 
missing  hydrogen  atoms  i.e.,  with  the  three  dangling  bonds  at  the  (111) 
surface. 

If  the  state  exists,  the  calculation  will  determine  its  location 
relative  to  the  top  of  the  valence  band,  that  is,  that  energy  required  to 
excite  an  electron  from  the  top  of  the  valence  band  into  the  surface  state. 

2.  INPUT  CALCULATIONS 

Certain  calculations  are  required  as  input  data  to  the  programs. 

These  include  the  tetrahedral  geometry,  coordinates  of  atoms,  symmetry 
point  group  analysis,  and  atomic  and  outer  sphere  radii. 

3.  THE  DIAMOND  STRUCTURE  (REFERENCE  31) 

The  space  lattice  of  diamond  is  face-centered  cubic.  A  primitive 
basis  of  two  identical  atoms  000;  1/4  1/4  1/4  is  associated  with  each 
lattice  point,  as  shown  in  Figure  3(a).  The  tetrahedral  bonding  of  the 
diamond  structure  is  represented  in  Figure  3(b).  Each  atom  has  four 
nearest  neighbors  and  12  second-nearest  neighbors.  Here  are  eight  atoms 
in  a  unit  cube.  The  diamond  structure  is  the  result  of  fourfold 
directional  covalent  bonding. 
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4.  The  C5H9  STRUCTURE 

In  the  C^Hg  structure  there  is  a  carbon  atom  at  the  center  of  the 
tetrahedron  formed  by  its  four  nearest  neighbors.  The  12  second-nearest 
neighbors  can  be  subdivided  into  four  groups  of  three  atoms  each,  each 
group  lying  in  a  plane  parallel  to  one  of  the  four  faces  of  the  tetrahedron. 
If  one  of  these  groups  is  removed,  the  corresponding  face  represents  the 
unrelaxed,  unreconstructed  (111)  diamond  surface.  The  remaining  nine 
second-nearest  neighbors  are  replaced  by  hydrogen  atoms  to  tie  up  the 
dangling  bonds,  resulting  in  a  C5Hg  structure  with  C3v  symmetry.  This 
structure  is  shown  in  Figure  4,  where  the  basal  plane  represents  the 
(111)  surface  and  the  dashed  lines  represent  the  three  dangling  surface 
bonds.  In  the  numbering  scheme  used  in  the  calculations,  the  carbon  atoms 
are  numbered  3,  4,  5,  8  and  9,  and  the  hydrogen  atoms  are  6,  7  and  10 
through  16.  "Atom"  1  denotes  an  outer  sphere,  concentric  with  atom  3, 
large  enough  to  encompass  the  cluster,  and  "atom"  2  denotes  a  missing 
H  atom  absorbed  on  the  (111)  surface  just  below  the  center  of  the  basal 
plane.  This  "vacancy"  atom  is  included  for  comparison  with  future  studies 
of  hydrogen  absorption  on  diamond.  With  respect  to  the  coordinate  axes 
shown  in  Figure  4,  the  coordinates  of  the  atomic  spheres  used  in  the 
calculation  is  given  in  Table  1,  both  in  units  of  aQ/4  and  in  atomic  units 
(au),  where  an  =  3.566  &  =  6.7387  au  is  the  length  of  the  cube  edge  in 
the  diamond  lattice  (at  0  K)  and  1  au  =  0.52918  A.  The  coordinates  of 
the  vacancy  are  determined  assuming  it  just  touches  the  center  carbon 
atom. 

5.  SYMMETRY  POINT-GROUP  ANALYSIS 

The  symmetry  point  group  of  C5H9  ,s  C3V  Rotation  about  the  basal 
plane  of  the  cluster  is  in  120°  steps,  such  that  the  identity  matrix  is 
generated  at  (C3)  . 

In  order  to  see  this,  the  rotation  of  axes  about  an  arbitrary  angle 
in  the  x-y  plane  is  shown  in  Figure  5. 
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Figure  5.  Rotation  of  Axes  in  the  x-y  Plane 


x'  =  x  cos  0  +  y  cos  (j  -  9) 

x’  *  x  cos  9  +  y  sin  0 

y'  *  x  cos(|  +  0)  +  y  cos  9 

y'  *  -x  sin  0  +  y  cos  0 


X' 

cos  9  sin  0 

V 

y'  * 

-sin  0  cos  9 

y. 
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2  3 

The  matrices  for  C^»  C^,  and  are  given  below.  Furthermore,  as 

can  be  seen  from  Figure  6,  triangular  symmetry  of  the  basal  planes  also 

evokes  reflections  depicted  by  the  a-axes.  Og  maps  position  seven  into 

position  six,  and  vice  versa.  Og  maps  position  five  into  position  seven 

and  vice  versa  which  maps  position  five  into  position  six  and  vice 

versa.  These  operations  generate  additional  symmetry  to  Cgy.  The  a 

operations  may  also  be  construed  as  permutations  and  illustrative  matrices 

2 

are  included.  Further  note  that  a g  =  OgCg,  that  is,  Og  can  be  generated 
by  rotating  the  system  240  degrees,  then  applying  Og  operation.  This  is 
similar  for  as  well  since  a ^  =  OgC^. 

6.  OTHER  INPUT  PARAMETERS 

The  radii  of  the  atomic  spheres  were  determined  by  assuming  touching 
spheres,  resulting  in  a  radius  for  atoms  3  through  16  of  one-half  the 
nearest-neighbor  distance  in  diamond,  or  1.4589773  au.  The  vacancy 
sphere  was  assigned  a  radius  equal  to  one-half  the  bond  distance  in  the 
ground  state  of  the  H 2  molecule,  or  0.7007067  au.  The  outer  sphere  radius 
was  taken  as  the  smallest  possible  for  encompassing  all  the  atoms,  or 
6.2239775  au. 

The  values  of  the  exchange  parameter,  a,  were  taken  as  ac  =  0.75928 
for  carbon  and  =  0.97804  for  hydrogen,  from  the  compilation  by  Schwarz 
(Reference  19).  For  the  inter-atomic  region,  a  was  determined  by  weighting 
a  and  «u  with  the  number  of  valence  electrons  (4  per  C  atom  and  1  per  H 
atom),  giving  aint  =  0.82661. 

The  maximum  value  of  t  in  the  partial -wave  expansion  of  the  symmetrized 
basis  functions  was  l  =  1  for  functions  constructed  from  the  hydrogen 
orbitals  and  i  -  2  for  functions  constructed  from  the  carbon  orbitals. 

The  fraction  of  the  new  potential  used  in  each  iteration  was  0.1. 

No  spin  polarization  was  assumed,  and  the  self-consistency  criterion  was 
a  1%  variation  in  the  potential  between  successive  iterations,  up  to  a 
maximum  of  50  iterations. 
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Figure  6.  Triangular  Axes  of  Symmetry  of  C^v 
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SUMMARY  OF  TWO-DIMENSIONAL  MATRIX  REPRESENTATION  OF  C, 


=  f1  °1 

0  ij 


-1/2  A/2 

C.  = 

J  -*''3/2  -1/2 


0^  “  0rCo  = 


-/f/2 

-  1/2 

a5  =  (o  -1. 

-1/2 

—✓3/2 

f-  1/2 

A/Z 

-✓3/2 

1/2. 

O7  *  OrC  —  = 

7  6  3  ^5/2 

1/2 

Trace  E  =  2 

2 

Trace  C3  *  Trace  C3  =  -1 
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SECTION  XII 

RESULTS  AND  CONCLUSIONS 

Application  of  the  theory  presented  in  Sections  I-IX  of  this  report, 
via  the  computational  scheme  outlined  in  Sections  X  and  XI,  resulted  in 
a  set  of  self-consistent  energy  levels  for  the  C^Hg  cluster  which  was 
convergent  to  at  least  one  part  per  thousand  after  50  iterations.  Table 
2  lists  the  occupied  and  low-lying  unoccupied  energy  levels,  and  their 
percent  change  A  between  the  49th  and  50th  iterations.  The  tabulation 
is  arranged  according  to  the  corresponding  irreducible  representations 
of  C^y,  namely  the  two  one-dimensional  representations  A1  and  A2  (maximum 
occupancy  of  two  electrons  in  each  level)  and  the  two-dimensional 
representation  E  (maximum  occupancy  of  four  electrons  per  level). 

Watkins  and  Messmer  (Reference  28)  made  a  calculation  similar  to 
ours,  but  on  the  cluster  CcH10.  They  were  interested  in  the  properties 
of  diamond  in  the  bulk  and,  to  the  extent  that  their  calculation  is 
representative  of  the  solid,  ours  is  representative  of  the  unrelaxed, 
unreconstructed  (111)  surface.  Comparison  of  the  two  calculations  is 
facilitated  by  making  use  of  the  correlation  table.  Table  3  (Reference  29). 

This  gives  the  splitting  of  the  levels  as  the  symmetry  is  lowered  from 
Td  for  C5H12  to  C3y  for  CgHg.  Thus  a  six-fold  degenerate  level  belonging 
to  the  T.j  representation,  for  example,  will  split  into  a  two-fold 
degenerate  level  belonging  to  and  a  four-fold  degenerate  level 
belonging  to  E. 

Table  4  compares  the  C5H12  valence  levels  from  Figure  2  of  I 

Reference  28  with  our  results  for  CgHg,  expressed  in  eV  (1  Ry  s  13.606 
eV).  The  comparison  Is  also  shown  on  an  energy  level  diagram  in  Figure  7. 

From  Figure  7,  it  Is  evident  that  the  cluster  corresponding  to  bulk 
diamond  has  a  well  defined  energy  gap  of  5.56  eV  between  occupied  and 
empty  states.  This  agrees  with  the  experimental  energy  gap  of  5.4  eV 
between  the  top  of  the  valence  band  and  the  bottom  of  the  conduction 
band  (Reference  12).  For  the  CgHg  cluster,  corresponding  to  the  presence 
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TABLE  2 


ONE-ELECTRON  ENERGY  LEVELS  OF  CgH9 


Level  Representation 

and  Occupancy) 

Level  Energy,  Ry 
(49th  Iteration) 

Level  Energy,  Ry 
(50th  Iteration) 

A(S) 

IS  core  of  C  atom 
#3.  (2) 

-20.523 

-20.522 

0.005 

IS  core  of  C  atom 
#4.  (2) 

-20.698 

-20.697 

0.005 

IS  core  of  C  atoms 
#5,  8  and  9.  (6) 

-20.675 

-20.676 

0.005 

A1  (2) 

-1.92808 

-1.92806 

0.001 

A1  (2) 

-1.57929 

-1.57857 

0.05 

A1  (2) 

-1.29193 

-1.29213 

0.02 

A1  (2) 

-1.18067 

-LI  8043 

0.02 

A1  (2) 

-0.85135 

-0.85174 

0.05 

A1  (0) 

-0.28917 

-0.28931 

0.05 

A1  (0) 

-0.15468 

-0.15409 

0.38 

A1  (0) 

-0.12237 

-0.12269 

0.26 

A2  (2) 

-0.96254 

-0.96374 

0.12 

E  (4) 

-1.52434 

-1.52489 

0.04 

E  (4) 

-1.14320 

-1.14330 

0.01 

E  (4) 

-1.02248 

-1.02345 

0.09 

E  (4) 

-1 .06616 

-1.06559 

0.06 

E  (1) 

-0.83327 

-0.83393 

0.08 

E  (0) 

-0.48203 

-0.48271 

0.14 

E  (0) 

-0.10953 

-0.10968 

0.15 

E  (0) 


-0.09989 


-0.09997 


0.08 


T~.  **:  •  ■  r.  r .  ■  ■  '  . 


AFVIAL-TR-83-1130 


TABLE  3 

CORRELATION  TABLE  FOR  Td  •*  C3V 


Td  reP* 


C3V  rep. 


A 

A 

E 

T 

T 


T 

2 

1 

2 


+  E 
+  E 


of  a  free  (111)  surface,  the  six-fold  degenerate  T ^  states  at  the  top  of 
the  valence  band  split  away  and  form  a  narrow  group  of  partially  occupied 
surface  states  of  width  0.24  eV  located  in  the  lower  half  of  the  gap. 

The  “Fermi  level"  in  our  cluster  model,  which  is  at  -11.346  eV,  lies  1.77 
eV  from  the  top  of  the  valence  band.  These  results  agree  quite  well  with 
self-consistent  pseudopotential  calculations  of  the  diamond  (111)  surface 
(Reference  8).  The  latter  obtained  a  band  of  partially  filled  surface 
states  approximately  0.2  eV  wide  and  1.7  eV  above  the  valence  band.  They 
also  deduced  a  value  of  the  work  function  <|>  of  about  7  eV  whereas  our 
model  predicted  4  ^  11  eV. 


In  summary,  the  SCF-Xa  scattered-wave  method  applied  to  diamond 
yields  results  for  the  electronic  structure  in  the  vicinity  of  the  band 
gap  in  very  good  agreement  with  available  experimental  data  for  the  band 
gap  in  the  crystal  and  with  band-structure  calculations  for  the  unrelaxed, 
unreconstructed  (111)  surface.  This  suggests  that  similar  cluster  models 
would  be  useful  in  studying  surface  properties  such  as  atomic  adsorption 
as  well  as  bulk  properties  such  as  point  defects,  in  other  covalently 
bonded  solids  having  the  diamond  structure,  e.g.,  Ge,  Si  and  GaAs. 


I 

i 
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TABLE  4 

COMPARISON  OF  VALENCE  LEVELS  FOR  C5H12  AND  CgHg 


CgHi2  levels  (eV)  and 
(representation) 

CgHg  levels  (eV)  and 
(representation) 

-0.56 

. 

(t2) 

-2.097 

(A-j ) 

- 

-2.22 

(A,) 

-3.936 

(A,) 

. 

a 

-6.568 

(E) 

|  -7.78 

(t2) 

-11.346 

(E) 

! 

-11.589 

(A,) 

-7.78 

<V 

-13.113 

(a2) 

-13.925 

(E) 

-9.44 

(E) 

-14.498 

(E) 

-10.56 

(t2) 

-15.556 

(E) 

-16.061 

(A,) 

-12.22 

(A,) 

-17.581 

(A-| ) 

-16.1 

(t2) 

-20.748 

(E) 

-21.478 

(A,) 

-20.6 

(A1 ) 

-26.233 

(A,) 
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APPENDIX  A 

sin(Kr  -  ^  +  N.) 

Proof  of  the  form  G^  ^ ^ - 

Solving  Equation  142  by  perturbation  theory: 

G;  <•  [K2  -  U(r)  -  G,  *  0 

1  r  (A-l) 


Letti ng 

=  g£  +  0  and  assume  that  the  product  JDU(r)  can  be  neglected, 
for  small  values  of  0,  our  perturbation  parameter.  (A-2) 


Substituting  Equation  2  into  Equation  1  we  have 

+  0"  +  [K2g0  +  K20  -  Ug0  -  U0  -  0]  =  O  (A-3) 


g”  +  [K2  -  U(r)  -  g^  +  0"  +  [K2  -  0  =  0 

r  r 

Looking  at  part  of  Equation  A-4  which  is  bounded  at  the  origin 


(A-4) 


+  [K2  -  ,  0 

dr  r 


g£  ^  sin(Kr  -  ^  £tt) 


(A-5) 


then  we  have 


q*<K2 

dr 


^7^-  )  0  s  Ufrjg^r) 


(A-6) 
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If  we  let  Equation  A-7 

*  ■  9,<r)  N(r)  (A-7) 

*'  -  Si"  *  SlN-  (A-8) 

»”  *  9p  *  29i«'  ♦  gtN"  (A-9) 

Substituting  Equations  A-7  and  A-9  into  Equation  A- 6  we  have 

9*N  +  2g'N'  +  g^N"  +  ( K^g^N  -  g£N)  =  U(r)9jl  (A'10) 

or 

N[g"  ♦  K2g£  -  Mttli  ♦  g(r  *  2g;N'  •  U(r)g£(r)  (A_n) 


Substituting  Equation  A-5  into  Equation  A-ll,  we  have 

g,N"  t  2g'N'  =  UMg^r)  (A.12) 

Multiplying  g^  into  Equation  A-12  we  have: 

( gf N ’ )  =  U(r)g|  (r)  (A-13) 

Integrating  Equation  A-13, 

g^N'  -  Jr  U(r)g^(r)  dr  (A-14) 

t+1 

Since  N'  must  be  bounded  at  r  =  0,  and  g£(r)  behaves  like  r  for 
small  r,  the  lower  limit  of  integration  must  be  zero.  Then 


For  large  r,  we  have  g£(r)  =  sin(Kr  -  y] 


(A-15) 
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or 


dN 

a?  * 


1 


[sin(Kr-  J0 


f  U(r)[g  (r)]wdr 
Jn  1 


s  esc2  (Kr  -  |°°  UCrKg^r^dr 


letting 


jl 


(r)j  dr  *  A. 


a  constant,  then 


3f  =  csc20<r  -  Az 


Integrating  Equation  17 

N  =  A^  ,/csc2(Kr  -  y-)  dr  =  -A£  cot(Kr  - 


since 


A  frntvi  -  d  /Cosx»  _  sinx  cos2x 
dx  Uotx)  dx  <117^  sinx  ~^T 


2  2 
-  sin  x  -  cos  x  _ 

- - — - -  esc  x 

sinx 


=>  /  esc  x  =  -  cotx 


then 


G*  =  gl  *  *  s  s*"(Kr  -  \  Itt)  -  [cot  ( Kr-y)  +  aU* 


x  sin  (Kr  - 
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(A-16) 


( A- 17) 


( A-l 8 ) 
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wbw?t&t*s  rvnwwvi&w  '/hw  nmwwjtvs*  pWiTi 


or 


.  a 

G£  =  sin(Kr  -  j  lit)  -  [cos(Kr  -  y-)  +  a  sin(Kr  -  -| 


otA,  A„ 


-  sin(Kr  -  J-  Air)  [1  -  -j£|  -  cos  (Kr  -  ~) 


letting 


1 


<*A  i 


s/> 


4>2  *[.-“> 


/V 

cosN£  and^sin  N£ 


X2 


aA 


*-.2 


s\/(- jf)  +  [1  -  — ]  sin  (Kr  -  ^  +  N£) 


cojp-  sin  (Kr  -  ^  +  N£) 


here 


sinN£  n.  N£ 


Ai 

X 


where  is  small. 
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APPENDIX  B 

PROPERTIES  OF  THE  FUNCTIONS  K^  (Kr)  AND  i  (Kr) 


I  •  T3‘  •  • (2&-1 


1  im 

K*1  ^ (Kr) 

r+O 

A 

1  im 

Kg1  ^ (Kr) 

r+O 

X. 

|/0 ) 
Vi 

II 

1. 

(Kr) 

(-1)1  **4^ 


dK^'lKr) 


forT"  '  JITT  + 

-  «<*>(»•)  *  K<»<Kr>  =  fiafe^ 
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\  (K«*) 


1-3-5 


••  ] 


i«(Kr) 

"TO" 


[Ai,  ,{Kr)  +  (A+l )(i.+1 (Kr)] 
2A+1  1  1  11 


\^r) 


dK0(1)(Kr) 


"d(KfT 


-  ^1}(Kr) 


diJl(Kr) 
d  (Kr) 


HI 

Kr 


A- 1 

T~ 
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APPENDIX  C 

THE  DERIVATION  OF  THE  INTEGRAL  SCHR0D1NGER  EQUATION 
The  Schrodlnger  Integral  equation  is  derived  as  follows: 

2m “  ^(r)  ’  6  ^(r)  +  u(^  *(r)  *  0  (C-l) 

or 

(72  +  K2)i{i  *  V4»  (C-2) 

where 

K2  a  2^E  and  V  =  2y  4- 
tr  n 

for  a  particular  solutfon  of  Equation  C-2  is  constructed  in  terms  of  the 
Green's  function  G(r,r‘)  which  is  the  solution  of  the  equation 

(V2+K2)G(r,  r')  -  -4tt6  (r  -  r'). 

The  expression  -  ^  /  G  (r,  r')  V(r')  ip(r ' )  d3r'  solves 

Equation  C-2  since  (V2  +  K?)<|>  *  Vty.  Multiplying  both  sides  by 
A  ^ 

G(r,  L' )  and  integrating:  /  dr*  G(V2  +  K2)^  *  /  V G  tf/dr'  since 
/  dr'  G V2t|>  *  /  i^VZG  dr',  /  dr'  |>V2G  +  K2G*]  =  /  VG  t//dr ’  = 

/  dr'  4>[V2G  +  K2G3 ■  -4tt/  dr*  rij<S(r-r')  =  ^(r)(-4ir) 

ip(r)  »  -  ~  /  V(r')G(r,  r')  ^(r1)  dr';  letting  * 

GQ(r,  r'),  then  ^(r)  »  /  V(r')  GQ(r,  r')  tl'(r')  dr'  as  was  to  be 
proved. 
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APPENDIX  D 

INTEGRAL  REPRESENTATION  OF  THE  FREE-SPACE  GREEN'S  FUNCTION 


Consider  applying  the  Fourier  transformation  to  the 


equation 


( v  +  e)  G(r  -  jr * )  =  -4ir  6(r  -  jr* ) 

i  K  JL  ( r  -  r 1 ) 

If  we  introduce  G{r  -  r/ )  =  /  g(K')  e  J  ~  ~  d3K ' 


iK'*(r  -  r') 

and  6(r  -  r')  *  — /  e  -  -  d3^  then  substituting 


these  expressions  we  have: 


/(V2  ♦  .)  g(K')  -  f>d3.=  -  I'W 


I  g(K')  [-K'2  ♦  e]  e1*'^ '  -Vic'  -  -V/elK''£'£d3  K' 

2tt 


Since  Gq(£ 


g(K‘ )  55  [-K'2  +  e]”1  •  hence 

2tt 


G(r  -  r* )  *  /  g(K')  eiK'’(r-r')d3K' 


—1  t  e^' ~  '  3 

■  -V/  - - p - dV 

2tt  [e-K'2] 

G(r  -  r ' ) 

r')  =  -7-5 - - ,  then 


1  ,  JK'-(r-r') 

GQ(r  -  r')  =  — Lj  /  ^ - 5 - 

0  '  (2tr)3  Ce  -  K'2] 


:,.xv.--.n^avvv,v^.v:.- yj 
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or 


(m7 


[e  -  K2] 


G0(r  -  L')  *  G0(r,  r',  e) 
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APPENDS  E 

DERIVATION  OF  AN  ALTERNATE  FORM  OF  GREEN'S  FUNCTION 
The  free-space  Green's  function  was  derived  in  Appendix  D.  It  is; 

G0(V  ^  f  d^P[ia-'  (S  *  ~  ^ 

0  a  S  (2ttV  J  e  -  q2 


where  R^  =  Rg  -  R*. 


For  the  case  where  8  =  ot,  the  alternative  partial  wave 


expansion  is  given  by 


i  exp(-K|r  -  r'|) 

^0(— *  r'i  e)  ■  -  TZ  - 

in-r*  I 


This  form  is  used  profusely  and  therefore  merits  a 
derivation. 

Let  us  consider  the  momentum  integral  (q  e  p)  representa¬ 
tion: 

Jg/r  2 

I  =  /  ~2 — ^  ^9-  w^ere  s  *  n  and  K  is  a  complex  number 
K  -q 


with  reK  >  0.  We  then  have 


108 
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■n  f 


2it  fit  eiqrcos0  2 


2 - g-  q  sin  9  d  9  dp  dq. 

K  -  q 

fit  iqrcos9 


rv  lqrcosw 

Looking  at  the  integral:  £-* — 7 —  sin  6  d  9  we  see 

'0  K-q 


Oiqrcos9 

that  it  is  just  =  -  - — * — ~- 
1qr(KZ-qZ) 


*  a  J.sinqr 
0  qr(K-qZ) 


K-q 


k  -  q 


where  -  ®  <  q  <  0.  The  Integral  -  f*  — 2"  ^9  the  principal 

K  -q 


2iri  / 

part  of  the  contour  integral  -  =y-  9f— ■ j  dq  in  the  complex  q 


plane. 


Here  the  integration  is  taken  over  the  Infinite  semicircle 
in  the  upper  half-q  plane  in  Figure  E-l . 


Figure  E-l.  Contour  Used  for  Evaluating  the  Green's  Function 
for  Case  (a) 
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The  contour  integral  is  equal  to  2-iii  times  the  sum  of  the 
residues  of  the  poles  of  the  integrand  which  lie  within  the  con¬ 
tour.  Hence  we  have 


Case(a)  im  K  >  0  -  pole  at  q  =  K 


27ria  for 


ffc 


iqr 


qTTK+q) 


«  2-r-  <K-» 


(«> 


-  2tt 


2  e 


iKr 


Therefore,  the  more  general  vector  representation 


i g;(r  -  r’) 


I  -  f  — 2 — 2 - d9.  becomes 

K  -q 


-  2tt 


2  eiK|r  -  n' 


I r  -  r' 


Similarly,  Case  (b). 


im  K  <  0  -  pole  at  q  =  -K  with  contour  of  Figure  E-2, 


r  pia-n 

gives  — 2  d9.  ■*'  -2it 

J  K  *q 


2  e 


-iKr 
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M 


Flqure  E-2.  Contour  Used  for  Evaluating  the  Green's  Function 
for  Case  (b) 


Therefore  the  more  general  vector  representation 


_  ,  o  0-iK|r  -  r1 

1  =  I  — 5 — 9 - d9.  becomes  -  2-rr  5 - 

K-qZ  (r-r'| 


Then 


.27r2  eiKl  rrt 
G0(ra’r6  je>  "  7TTT 


(2ir)  |r-r'| 


elKl-  "  H' 

4n 


r  -  r’ 


as  was  to  be  proven. 


Ill 
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APPENDIX  F 

DERIVATION  OF  ASYMPTOTIC  FORMS  FOR  HANKEL  AND  BESSEL 
FUNCTIONS  OF  REAL  ARGUMENT 

Consider  the  equation  (V2  +  KgMr)  =  S(r),  where  V2  is  the 
spherical  coordinates  Laplacian. 

A  solution  exists  when  one  finds  the  Green's  function 


G (r,  r')  -  (V2  +  K^)'1  5(r 

and  then  uses 

Mr)  3  / G(r,  r')  S(r')  dr' 


-  r') 


(F-I) 


(F-2) 


Since  is  complex  in  general  and  the  eigenvalues  of  VZ  are  real, 
the  operator  (V2  +  Kq)  cannot  have  a  zero  eigenvalue.  Hence. the 
existence  of  the  Green's  function  G  is  a  possibility. 

In  order  to  solve  this  expression  and  determine  the  form  for 
the  Green's  function,  some  preliminary  results  have  to  be  obtained 
We  know  that  6(r  -  r')/3  {^-)3  /e1-  ^  <*K  that  is: 


6(x-x' )  <5(y-y')  fi(z-z') 


(F-3) 


-4  rri 

(  2tT  )  J  -00  J  -00  j 


-H»  f4»  ,*»  i[K  (x-x* )  +  K  (y-y‘ )  +  K  (z-z’) 


dKxdKydKz 
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Transforming  to  spherical  coordinates,  we  have  (see  Appendix  G) 


6(r  -  r')  =  g(r-r,)6(fl-Q,)6(»-»<) 
rZsin9 

Thus 


r  r 

r  sln  (2tt)3  '0  JO  Jo 

i K* (r-r* )  2 


(F-4) 


(2tt)' 

K“  sin  adK  dadB 


The  expansion  of  a  plane  wave  e^— — ,  traveling  in  an  arbitrary 

direction  a,B  (i.e.  K  has  components  K  *  K  sinacosB;  K  = 

«  y 

K  sinctsinB,  Kz  *  K  cosa)  is 


iK-r 


+fc 


•  l  l  i  Oo(Kr)  Yj  (a.B)  ¥"(0,0) 

ll30  m=-t  *  * 


(F-5) 


Substituting  Equation  F-5  into  Equation  F-4  for  e1-  -  ancl  tlle  comPlex 

j  1/  ,  ^  • 

conjugate  of  Equation  F-5  for  e  into  Equation  F-4,  we  have: 


S(M., 


•I  ,  r2l T  rTT  roo  co  +1 

*  L  L  L4”  r  i  <“-s) 

J0  J0  J0  1=0  m=-i  *  * 


(F-6) 


oo  +]. 


X  yJ(9,0)  •  4tt  l  l  (-i  )LjL(Kr’  )Y^(9,0)Y^(0  '  ,0 ' ) 
1*0  M=-L  L  L  L 


x  K  sina  da  dB  dK 
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or  in  vector  form 

Mr-r')"|  I  l  Hl(-1)L  vf  (9',*1)  Y?  (9,8)  .. 

11  £,m  L,M 

x  |  j)l(Kr)  jL  (K r ' )  K2dK  j  j  Y™(a,8)Y^(a,8)  sina  da  dB- 
Because  of  the  orthonormality  condition: 

j  J  Yj  (a,B)  y[*  (a,B)  sina  da  dB  =  6£L  <5^  (F-8) 

Therefore, 

6(r  -  r')  =  |  l  Y?  (0*,4>*)  Y?  (0.4>)  f  j£(Kr)  j  (Kr')K2dK 

Z,m  >0  (F-9) 

2  2-1 

looking  back  at  Equation  F-l ,  namely:  G(r,  r,')  =  (V  +  KQ) 

<5(r  -  r' )  we  see  that  the  Green's  function  will  have  an  expansion 

in  terms  of  spherical  harmonics  and  Bessel  functions.  Therefore 

we  can  let 

G(r,  r')  =  l  g£(r,  r ' )  yJ  (Q',*'}  Ym  (0,$)  (F-l.O) 

Z ,  m 

It  is  necessary  to  address  the  problem  of  determining  the  recip- 

2  2 

rocal  of  the  operator  V  +  Kq.  To  do  this  let  us  consider  the 
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2  2 

operator  -a  V  +  1  in  spherical  coordinates,  and  the  inhomogeneous 
partial  differential  equation, 


Lf(r)  =  (V2  +  Kq)  f(r)  =  s(r) 


(F-ll ) 


where  r  denotes  the  vector  with  components  r  sin<|>  cos  e,  r  si n<f>  sin<j> 
and  r  cos<t>  and 


v22^irr2^.+  1 


r  3r  3r  r2sin0 


3  sinQ  4  +  -J-  g—  X 
36  00  r2sin20  *** 


(F-12) 


Now,  g(r,  jr ' )  =  lT1  <5(r  -  r ' ), 
gives 


G (r,  r')  *  L"1  f  I  Yj  <0V<J>')  yJ  (0,<j>)  j~  ^(KrJj^Kr *  )K2dK 


Jt,m 

__  m  r  j»(Kr)  .1  p ( Kr ' )  z 

l  Y?  ;0,,<J>1)  Y?  (0  ,d>)  1  - K - K^dK 

£,m 


(F-13) 


i — r 

0  Kq  -  & 


where  the  form  of  the  denominator  can  be  seen  from  Appendix  D. 


Comparing  Equations  F-10  and  F-13  we  have 


g^r,  r') 


2f 

*  Jn 


J^Kr^fKr')^ 
0  K I  -  K2 


(F-14) 


Since  g^(r,  r')  is  even  function  in  K,  the  integral  can  be  re¬ 
written  as 
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g0(r,  r') 


!  f  00jJl(Kr)jJl(Kr * )  ICdK 


H 


1 — r 
Ko  -  K 


(F-l 5) 


where  the  limits  have  been  changed  from  0  <  K  <  °°  to-~  <  K  <  °° 
We  can  write  in  partial  fractions: 


B 


K0-K*  ’  Ko  *  K  K0  -  K 


Thus 


K2  =  A(Kq  -  K)  +  B(Kq  +  K).  For  KQ  =  K,  K2  =  2Bk,  and 


B  =  For  Kq  =  -K,  K2  =  -2KA  and  A  »  -  £  . ' 


K 

“2 - 2 

K0  "  K 


Now 


-K/2  +  _K/2 


Kq+K  Kq-K 


d< 

2  K-KQ  '  K+Kq 


1  +  1 


K  +  K 


2K 


K  +  KQ  '  K  -  Kq  "  K  -  Kq  +  +  Kq  "  K  -  KQ‘ 


(F-l  6) 


Since  the  term  in  parentheses  on  the  right  is  an  odd  function 
in  K,  it  will  make  no  contribution  to  the  Integral.  Therefore, 

hlr-  r'> =  -  ?  f  h^hiKr">  K-nc dK 

*  -oo  n 


(F-l 7) 


( F-l 8 ) 


t  -  1  r 


We  shall  integrate  Equation  F-l 8  by  the  use  of  contour  integration. 
Since  j'A(x)  *  \  [h^(x)  +  h|2^(x)]  and  hj^(x)  ^  ~  . 

(2)  e"ix 

h^  ;(x)  “v  — —  for  x  *  ,  the  integral  will  consist  of  terms 


116 


'~1  -1  4.C4JttLl 


AFWAL-TR-83-1 1 30 


behaving  like  e*^r  and  e*^r  .  For  the  case  r  >  r' 


/  \  _  1  f*  1  iKr  -iKrw  iKr'  -iKrS  v 

9t(r,r  )  -  -?J_?  (e_LS - ][fi - ±* - ](A  ) 


X  r  fJKr .  oKr*  v 

2it  J— l  Kr2  )  K^cV0"'* 


Then  the  two  parts  of  the  term  jA(Kr)  determine  how  the  contours  are 


to  be  closed.  The  h^(Kr)  part  of  ,i ^ ( Kr )  gives  a  contribution  at 


K  -  Kq  when  the  contour  is  closed  around  the  upper  half-plane. 


Since  the  h^'(Kr)  part  causes  the  contour  to  be  closed  around  the 


lower  half  plane,  where  no  pole  is  included,  it  gives  zero  contri¬ 


bution. 


•  g»(r,  r')  =  =fj|t  W1)  h^'1  (KQr)  for  r  >  r 


-i'W')  h^1  }(K0r) 


A  similar  calculation  in  the  case  where  r  <  r'  yields 


9jl(r,  r')  =  -iKQ  hju  (K0r* )ja  (KQr)  r  <  r' 


We  now  give  the  details  of  the  contour  integration.  We  can  write 


i  n\  io\  i  iKr  -iKr 

Jl<*r)  ■  7  (h‘,,(ltr)*h<2)(Kr))  -  \  <V  +  V~> 


2  v  Kr  Kr 


v  '  •  v  *  ’ 
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-1  r  e^Krj.(Kr')dK 
•'• 9  (r>  r,)  — 


±C 


e^j^Kr'JdK 


«  r(K-K- 


Call  these  two  terms  -1^  and  -I2  respectively:  nA(r,r')  =  -I,  -  I2- 


1  f  e  r 

Looking  at  integral  I,,  we  seek  the  solution  to  -  0  — — - 

'c  0 

1  I  eUi  (c,z)d2 

where  Z  *  Kr;  Zn  -  Knr;  let  — *ot.  Then  we  have  —  0  _ 

u  u  r  r  Jr.  z-7. 


Consider  the  contour  shown  in  Figure  F-l. 


Figure  F-l.  Contour  for  Evaluation  of  g  (r,r') 


ruwvwuLV't.’vj’w^Bw.!  m  v.'~.,..rr  -a».,  g-.- 
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Integration  around  this  contour  gives 
iz.(az)dz 

"•  ItV-I 


iz.(az)dz  ix.(ax)dx 

6  h  V6  e  ^ 


X  -  X, 


0  J  -R 

+  |0  h1^*0*6  \4(a  [xo+eei§])ieei0d0 
*n  (xQ+ee^)  -  xQ 

/R  i^x  +e) 

+  (  *■  ji(«[x0+e])dx 

!*0+*  (x+e)  I  xQ 

iRe10  ,Q 

f-R  e  jA[aRe10]iRg0d0 


Re”  -  2, 


Taking  the  limits  as  R  -*•  «  and  e  ■*  0,  we  see  that  the  first  and 

third  integrals  combine  to  give  the  integral  [  e  jftfcxjd* 

x  -  x. 


The  second  integral  is: 


lim 

e-*0 


e^(xn“®®  )  iffl  \A) 

jo  Ca(xn+ee10)]ieei^ajD 


(x0  +  ee1®)  -  x0 


pi  xn 

*  **ine  0  j£(axQ) 


and  the  fourth  integral  is 


\ 
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11m  f 
R-*»  J-f 


lim  f 
R-*°°  J-f 


1R410 

*  j^CoiRe10]  1Rei0dQ 


Re1 6 


-  z. 


ie1Rcos0  e-Rs1nQ  j  dQ  ^  Q 


since  j.[aRJl’0]  Is  a  bounded  function. 


I 


e1Zj0  C“z]dz 


r- 

J  -oo 


z  -  z. 


^iXjA[ax]  dx 


1x0 

0  *  -Trie'  0  j£(otx0) 


x  -  xo 


eixj£[ax]  dx 
x  -  x„ 


ixo 

■Hrie  0  j£(ax0) 


1 K  f* 

changing  variables  back  from  z  to  K  gives  +  te  0  J £, ( K0r * )  for  the 
value  of  this  integral. 


For  I2,  we  can  change  sign,  letting  -Kr  =  z  to  give  the  inte- 
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2irr 


1 

2vr 


k 


r 


e,xJt(-*r')d* 


■X  -  X, 


e<xJt(-*r')d, 

X  +  XA 


Since  the  pole  x  -  -Xg  does  not  exist  on  the  aforementioned 
contour,  then  integral  1^  ■  0. 


'i 


eixja[dx]dx 


0 


(1), 


00  x  -  x. 


=  -Trie  uJA(ctxQ)  =  -irih^  '(x^j^axg) 
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APPENOIX  G 

COORDINATE  TRANSFORMATION  OF  DIRAC'S  DELTA  FUNCTIONS 

A  locally  integrable  function  f(x)  *  f(xj,...,xn)  defines  a 

distribution  through  <f,0>  *  /R  f(x)  0(x)  dx. 

n 

It  Is  useful  to  change  coordinate  systems  at  times.  To  show 
how  this  Is  done,  we  let  Uj,...,Un  be  new  coordinates  that  can  be 
obtained  from  x^,...,xn  by  the  transformation  law  U  *  U(x),  that 
is,  U.|(x],...,xn),...,Un  *  Un(x^,...,xn).  If  the  Jacobian  J  of  x 
with  respect  to  u  is  positive  everywhere 

> 

<f,0>  =  /  f (x)  0(x)  dx  (G-l ) 

Rn 

=  /  f (x  (U) )  0(x(u))  Jdu  =  J  f  (u)J0(u)du 

u  space  u-space 

**» 

where  0  =  0(x(u));  f(u)  =  f(x(u)). 

Equation  G-l  can  be  used  to  Interpret  coordinate  changes  for 

arbitrary  distributions.  Considering  f  *  6 (x-x * )  =  6(x] -x ' ) _ 6(x  -x'  ), 

then  <f,4»  =  ♦(x ' )  =  4>(u'),  where  u'  =  u(x').  Thus,  for  Elation  G-l  to 
be  consistent,  we  need 

f (u)J  =  dfUj-Uj)  ...  <S(Un-U^) 

or  since 

f(u)  *  5[x(u)-x'] 

we  f 1 nd 

6(UrUj)  ...  6(Un-U’)  (G_2) 


6(x-x ' )  = 


J 


v  v <m  -•■-.■  ^  v  tot  ^ 
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If  we  only  know  that  J  ^  0,  Equation  G-2  can  still  be  use,  that  is, 
if  J  does  not  vanish  at  U‘,  because  only  the  value  of  J  at  u'  affects 
the  right  side  of  Equation  G-2. 

To  calculate  the  transformation  of  the  delta  function  from  x-space 
to  spherical  space,  we  are  required  to  calculate  the  Jacobian  for 


cose. 


.  3(x,y,z)  _ 
a(rt0,0 


«(♦- 

<  =  r 

ax 

ax 

ar 

50* 

nr 

& 

3* 

ar 

as 

a? 

az 

az 

az 

ar 

30 

nr 

sinQ  cos0  r  cosQ  cos0  -  rsinQ  sin0 
sinQ  sin0  r  cosQ  sin0  rsinQ  cos 0 
cosQ  -  r  sinQ  0 


*  r2cos20  sinQ  cos20  +  r2  sin3Q  sin20 
+  r2  sinQ  co$20  sin20  +  r2  sin30  cos20 
=  r2  cos20  sinQ  +  r-  sin3Q  =  r2 


sinQ 
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APPENDIX  H 

SECULAR  DETERMINANT  AND  EIGENVALUES 

To  find  energy  eigenvalues  for  a  state,  the  matrix  elements  are 
constructed  and  the  determinant  is  examined  as  a  function  of  the  energy 
to  find  the  zeros.  A  typical  plot  of  the  secular  determinant  for 
crystalline  zirconium  is  shown  in  Figure  H-l .  The  roots  are  found  either 
graphically  or  by  interpolation  schemes. 

The  determinant  of  a  large  matrix  is  evaluated  by  a  sweep-out 
procedure  which  results  in  a  matrix  with  only  zeros  on  one  side  of  the 
main  diagonal.  The  determinant  for  this  kind  of  matrix  is  simply  the 
product  of  the  diagonal  elements.  The  procedure,  used  in  our  energy 
program,  can  best  be  illustrated  by  an  example. 

Given  the  matrix 

M1!  = 

we  add  (-2)  times  each  of  the  elements  of  row  1  to  the  corresponding 
elements  in  row  2.  Similarly  we  add  (-3)  times  row  1  to  the  third  row. 

The  result  is 

'2  4  6' 

0-2  4 
0  4  14 

* 

which  has  zeros  under  the  first  diagonal  element.  The  same  procedure  is 
applied  to  the  second  diagonal  element.  We  add  (+2)  times  the  elements 
in  row  3  to  get 

'2  4  6 
0-2  4 
0  0  22 


2  4  6 
4  6  16 
6  16  32 


the  determinant  (2)  (-2)  (22)  *»88. 
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The  operations  on  a  particular  matrix  element  which  result  from 
applying  this  sweep-out  procedure  to  the  entire  matrix  are  given  by 

m1u  /  muu 

U*1 

where  I  is  the  smaller  of  1  or  j.  M1J  is  the  original  matrix  element  and 
m^  is  the  corresponding  matrix  element  after  the  sweep-out  procedure. 


0.1 


<S1 

o 


c 

s 


S  -0.1 

i. 

<o 

'a 

3 


-0.2 


Figure  H-l.  Secular  Determinant  as  a  Function  of  Energy  for 
Zirconium  at  the  Symmetry  Point  r 


Energy  (Ry) 
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